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Abstract 

This project developed fluid circulation models for the two- and three-dimensional 
Lagrangian shallow water equations. There were two stages to this development: in the 
first, the two-dimensional shallow water equations were transformed from first principles 
of oceanography into a serial implementation in MATLAB. In the second part, a serial 
implementation of the three-dimensional shallow water equations, developed by Dr. James 
Greenberg, was modified to run in parallel on many nodes of a computing cluster. 

The serial, one-dimensional model includes one lateral degree of freedom—the re¬ 
direction. The vertical or z-direction is modeled by layers; velocity in this direction was 
removed by a series of transformations to the governing equations. Development started 
with the conservation of mass and conservation of momentum equations. The traction terms 
in these equations were approximated using a method previously established by Mellor and 
Blumberg in their work on the Princeton Ocean Model (POM). These equations were scaled 
and then transformed to an s-coordinate system, again following the established method 
of POM. Once a set of partial differential equations were derived, these equations were 
discretized in space and time and solved in MATLAB. The implementation of the model in 
MATLAB allows the user a wide range of initial conditions for factors such as the bathymetry, 
initial area covered in fluid, magnitude of friction coefficients, and more. For a given set of 
initial conditions, the model steps forward in time by user-designated time steps solving 
for position, velocity, and depth of the fluid at each step and visually representing this 
information with appropriate graphs. 

The next stage of the project jumped forward to a serial three-dimensional implemen¬ 
tation of the shallow water equations developed by Dr. Greenberg. This model is similar to 
the two-dimensional model but with the added complexity of two lateral degrees of freedom— 
both the x- and y-directions—while the vertical component is still treated the same as in 
the two-dimensional model. In order to split the two-dimensional domain grid into a set of 
smaller domains assigned to the various processors, the special geometry of that grid had 
to be taken into account. Once a scheme was in place to divide the grid while maintaining 
all of its special properties, the computation on each sub-domain was performed with the 
same program that operates on the entire domain. This allowed easy implementation of 
a parallel solution: each node ran a modified serial implementation on a subsection of the 
larger problem that had been carefully separated from the whole. A process was created to 
allow the nodes to communicate data from the edges of their sub-domains as they advanced 
forward through time. 

There are two deliverable products from this project. First, there is a serial two- 
dimensional model of fluid circulation that takes into account many different user-designated 
initial conditions and can be useful for determining how well the mathematics of this model 
can approximate physical phenomena. Secondly, this project produced a three-dimensional 
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parallel model that serves as a proof of concept for future development of more advanced 
parallel models. 

Key Words: shallow water equations, free boundary, Lagrangian, s-Transformation, 
coastal ocean circulation model, MATLAB 
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List of Variables 

x, y , z: Directions in the Cartesian coordinate system. 

X(t), Y(t ), Z(t): Cartesian components of particle location in Lagrangian representation. 

s: Coordinate to which the Cartesian z is transformed, s is 0 on the bottom surface of a 
fluid and 1 on the upper surface. 

u(x,y, z,t): Physical velocity in the x direction in Eulerian representation. 
v(x,y, z,t): Physical velocity in the y direction in Eulerian representation. 
w(x,y, z,t): Physical velocity in the z direction in Eulerian representation. 
p : Density of the fluid. Assumed to be a constant equal to 1. 
p(x,y, z,t): Pressure of the fluid. 

/: Coriolis force. 
g: Acceleration due to gravity. 

a a p. Traction on an elemental cube, the force is acting on face of the cube a in direction b. 
The numbers 1, 2, and 3, correspond to directions x, y, and z. 

E: Eddy viscocity. 

u, v, w, p: Non-dimensionalized forms of u, v , w, and p; each is a function of x, y, z, and t 
H : Vertical distance scaling factor. 

L : Horizontal distance scaling factor. 

U : Horizontal velocity scaling factor. 

A: Aspect ratio; defined by 

e: Inverse Reynolds number; defined by jjp. 

C 2 : Inverse Froude number; defined by r j^. 

Kf r : Coefficient of bottom surface friction. 

K w : Coefficient of upper surface friction; namely, wind friction. 
u\ Horizontal velocity when transformed to s-coordinates. 

0: A mathematical expression used to simplify the governing equations. Has the value 
0 = w — ( a x + sh x ) u — sh t . Can be used in post-processing to determine values for w. 
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u: Depth-averaged component of horizontal velocity; depth indepedent. 
u: Residual component of horizontal velocity; depth dependent, 
u: Vector of horizontal velocities at all marker points in the fluid domain. 
u w : Vector of wind velocities at the surface at all fluid marker points. 

M: Mass matrix. 

M: Viscosity matrix. 

M w , Mf r : Upper and lower surface frictions matrices. 



1 Introduction 


The United States Navy has enjoyed years of dominance in the open ocean. It is 
unchallenged as a blue-water power, but is facing increasing operational challenges in the 
littoral regions of the world. The majority of the world’s population resides in coastal regions, 
which id also the scene of most of the world’s conflict areas. Therefore, it is important for 
the Navy to expand its operational capability in shallow water regions such as estuaries and 
enclosed bays. 

One significant way in which littoral regions differ from the open ocean is in the 
motion of the water there. In the open ocean, currents are largely predictable.The different 
seasons’ prevailing currents have been charted since the days of the earliest oceanic explorers. 
However, currents in coastal regions are heavily affected by the shape of the local terrain— 
known as bathymetry—as well as tides. While in the open ocean the sea floor is of little 
consequence to all but the deepest-diving submarines, in coastal regions the sea floor affects 
ships of all sizes. Tidal influences in littoral regions affect navigation for both surface ships 
and submarines, as well as amphibious landings, special warfare, and non-combat operations 
such as preparation for storm surges. 

To deal with these challenges more appropriately, the Navy would like to be able to 
simulate the motion of the sea in shallow water regions. To enable this the Navy needs a 
model that can predict the state of the water, or fluid, in the future based on the bathymetry 
and conditions in the fluid at the current time. The inputs to such a model would be the 
initial conditions of the state variables at the current time. The state variables are the 
unknown dependent variables of the problem; quantities such as velocity of fluid particles, 
height from the lower to the upper surface of the fluid, or water column, pressure, density, 
and more. The independent variables of this problem are location in three dimensional space, 
commonly expressed in Cartesian coordinates by the triple ( x , y , z), and time. The output of 
this program will be the predicted values of these state variables at some time in the future. 

This model is created using standard oceanographic equations as the starting point 
and develops these into the governing equations. The system of governing equations has two 
sources of non-linearity: the bathymetry and the Eulerian form of the governing equations. 
Because the problem is non-linear, a solution must be found through numerical methods. 
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2 The Princeton Ocean Model 

Computer models already exist that predict future values of state variables from 
information available at the current time. One important model is the Princeton Ocean 
Model (POM). Much of the Navy’s coastal ocean modeling software is derived from POM. 
The Princeton Ocean Model is a pioneering ocean modeling software developed by Alan 
Blumberg and George Mellor in 1977. POM is a numerical model that is not grid specific-it 
can be applied to any bathymetry. It is written in FORTRAN and has been extended and 
modified significantly over the past thirty years. The mathematical foundation of POM 
is a turbulence closure model developed by Mellor in 1973 and expanded in collaboration 
with Tetsuji Yamada over the next ten years. The Mellor-Yamada model is based on older 
turbulence closure hypotheses developed by Rotta [11] and Kolmogorov [5]. The governing 
equations of POM evaluate fluid properties, including velocity, temperature, and salinity, in 
a three dimensional space corresponding with Cartesian rectangular coordinates [10]. 

POM is developed from first principles of oceanography: conservation of mass, con¬ 
servation of momentum, heat transport, and salt transport. These equations resemble the 
Navier-Stokes equations, except that they are relevant to turbulent rather than laminar flows. 
For shallow water estuary problems, all flow is assumed to be turbulent. POM replaces the 
Navier-Stokes stresses that are applicable to laminar flows with Reynolds stresses that are 
used for turbulent flow. 

The model we propose in this project uses a similar mathematical approach to POM’s 
leap-frog system for computing mean velocities. Both POM and our model first update the 
mean velocity at each point based on data from previous time steps, then update particle 
positions using the just calculated values of the mean velocity. After these steps, other 
quantities can be computed and then the cycle repeats by updating mean velocity at the 
next time step. This approach helps to preserve conserved quantities. 

POM was not intended to address shallow water and estuary problems; rather, it was 
originally designed to solve for fluid motion on an oceanic scale and was later modified to 
deal with shallow water regions. It operates on horizontal scales of between 1 and 100 km, 
while the vertical scale is in tens to hundreds of meters. The time scale of POM ranges from 
tidal to monthly intervals [8]. 

2.1 Eulerian and Lagrangian Representations 

The Eulerian and Lagrangian approaches are two different representations for fluid 
motion. In the Eulerian or field representation, the velocity and acceleration of a particle are 
represented in terms of that particle’s position in space. In the Lagrangian representation, 
the velocity and acceleration are associated with specific particles. The Lagrangian is called a 
particle tracking representation because each individual particle has velocity and acceleration 
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functions assigned to it that follow that particle through all positions and times. In the 
Eulerian representation, velocity and acceleration are functions of the position in space at 
any given time [7]. 

Both POM and our model develop from the first principles of oceanography repre¬ 
sented by Eulerian field equations because these field equations are simpler to develop than 
their Lagrangian counterparts. POM remains a primarily Eulerian model even in implemen¬ 
tation. Since it is based on held equations, it does not do well with certain free boundary 
problems because the original model was not designed to track these boundaries. The free 
boundary is the physical interface between land and water. For a wet-dry problem, the model 
is asked to find the times when a surface is either wetted or dried because of tidal action, 
a storm surge, or other physical phenomena. The Eulerian field equations do not explicity 
track the motion of particles from a permanently wet area to an area that is sometimes dry. 
Instead, POM tries to calculate the velocity of the fluid at fixed grid points regardless of 
whether those points are in the wet or dry area. Furthermore, it has no mechanism to mark 
a “wet” area that has become “dry” or vice-versa [1]. 

Since POM has trouble distinguishing between wet and dry areas, and because we are 
interested in solving this sort of free-boundary problem, we take a slightly different approach. 
We start with the same equations, and develop them in a similar Eulerian method. However, 
we create a Lagrangian grid covering the lateral space of the fluid domain that moves with 
the flow of the water. The lateral space of the domain are the dimensions other than the 
vertical dimension. For example, when implementing a two-dimensional model from the two- 
dimensional shallow water equations, the lateral space is along the ^-coordinate while the 
^-coordinate is the vertical space. With the three-dimensional model, however, the entire 
x — y plane is the lateral space. The nodes in our grid are called fluid-markers and they 
are not fixed in space, nor are they fixed relative to other nearby fluid-markers. Instead, 
each fluid-marker point is associated with a particle and markers move as the velocity of the 
corresponding particle is re-calculated at each time step. 

This method is better able to track the wet-dry boundary. By defining a certain set 
of fluid markers as the “edge” of the water, we can track where these particles move and thus 
track where the free boundary moves, allowing us to easily distinguish between wet and dry 
surfaces. We believe that for certain types of problems, especially those that deal with the 
shallow water estuaries in which we are interested, this free boundary tracking is a better 
way of accurately modeling coastal ocean circulation. 

2.2 ^-Coordinate Transformation 

The most important development of the Mcllor model is the introduction of the a- 
coordinate system which helps in dealing with significant topographical variability. The a- 
transformation replaces the depth coordinate with a, a number that ranges between —1 and 
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0. Each value of a has a different ^-coordinate at different points (x, y). This transformation 
reduces the vertical component of the domain of the problem to a regular rectangular shape 
in <7, and pushes the complexity added by irregular bathymetry into the governing equations. 
Since the depth dimension is divided into many “layers,” each corresponding to a numerical 
value of a, it is fundamentally different from the x and y directions [1]. 

The model we are developing uses a variation on the ^-coordinate transformation 
that we refer to as an s-coordinate transformation. Where a ranges between 0 at the upper 
surface of the water and —1 at the lower surface, s varies from 1 at the upper surface to 0 
at the bottom surface. 

2.3 Numerical Viscosity 

While POM has a complicated mechanism that uses the Reynolds Stress terms from 
the physical equations to remove high frequency oscillations from the model, we have chosen 
a simpler method. We approximate POM’s eddy viscosity (e) with a numerical viscosity. 
Because we make this approximation, we have the freedom to set the value of numerical 
viscosity and can therefore choose a value that will help cancel terms in our governing 
equations. We nse a value of e that is dependent on both the size of the time step and the 
size of the interval between marker points to help our numerical model to converge. 

2.4 Asymptotic Solution to the Velocity Residual 

POM associates a location in space represented by the point ( x , y) with the averaged 
properties of the infinitesimally thin column of fluid above that point. The velocity at all 
depths above (x, y) is averaged into a mean velocity associated with that point. The dif¬ 
ference from the mean at each vertical position in the water column is called the velocity 
residual. POM attempts to directly solve the vertical profile equation for this residual. How¬ 
ever, we have made a significant contribution to work in this held by creating an asymptotic 
solution for the vertical profile. We show that the solution to the vertical profile can be 
split into transient and steady-state components and by obtaining appropriate estimates, 
demonstrate that the transient component can be ignored because it rapidly decays to zero. 
Therefore, the vertical profile can be quickly approximated with our model. 
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3 Physical Equations 

The derivation of equations in this and the following sections borrows from many 
sources; including [2], [4], [6], [7], [9], [12], and [13]. 

Both our model and POM start with the basic equations of physical oceanography. 
The scale of the shallow water regions for which we design this model is such that the /3-plane 
approximation can be adopted. Using the /3-plane approximation, we model the curvature 
of the earth with a linear approximation and a Cartesian coordinate system; this simplifies 
the mathematics of our model. The cost of this assumption is that a linear factor must be 
added to account for the motion of the spherical earth. This is the Coriolis effect and will 
be discussed in greater detail later in this section. 

We use both held and particle tracking equations as we develop and evaluate the gov¬ 
erning equations of this model. When a particle is tracked in a Lagrangian representation, its 
location in space is described by ( X(t ), Y(t), Z(t )) while the velocity components associated 
with a point in space in the Eulerian representation is (x(t), y(t), The velocity of a 

particle in either representation is described by the vector V(x, y, z, t ) where 

V = i u + j'f; + ku; 

is such that u, v, and w represent the component of the velocity in the x, y, and z directions, 
respectively. The identities for u, v, and w in Lagrangian notation are 

^X(i)=u(X(i),Y(t),Z(t),i), 

at 

A'(i) = v(X(i),Y(i),Z(i),i), (1) 

at 

^Z(f)=w(X(i),Y(t),Z(i),i). 

at 

Using these notations, we derive our governing equations. In (Figure 1) we consider 
an elemental cube with its faces aligned to the coordinates (x,y,z). If this cube is small 
enough that u, v, and w can be considered constant on all its faces, then the net volume 
how through the x and x + Ax faces is 

(u + UxAx) AyAz — uAyAz = UxAxAyAz. 

We use a similar derivation for the y and z faces and sum the net volume how through all 
six faces to obtain 

(ux + Vy + Wz) Ax AyAz. 

The Boussinesq Approximation states that density diherences between fluid elements can be 
ignored when approximating a how, except where the density appears in a term multiplied 
by g, the acceleration due to gravity. This approximation applies in the case of the equation 
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above. Since we assume that fluid density is constant, we can assume that the mass inside 
the constant volume AxAyAz does not change. Therefore, we see that 

Ux + % + wz = 0 ( 2 ) 

or in vector notation 

V- V = 0. 

This is the equation for mass continuity, or conservation of mass, which constitutes one of 
the four governing equations of our model. 

The other three equations derive from Newton’s Second Law of Motion. These mo¬ 
mentum equations balance forces with accelerations. For the elemental cube in (Figure 2), 
the net force in the x-direction caused by pressure from the adjacent fluid is 

—pxAxAyAz 


while the net force due to viscous or turbulent stresses in that direction is 


((dn)x + (d -21 )y + (<7 3 i)i) AxAyAz. 

The first subscript on the stress symbol a signifies the coordinate normal to the face of the 
cube on which the stress acts, and the second subscript is the direction of the stress, where 
we replace x with 1, y with 2, and 5 with 3 to avoid confusion with the partial derivative 
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Figure 2: An elemental cube showing force in the ^-direction 


subscripts. These two forces are equated to the product of mass and acceleration according 
to Newton’s Second Law. The mass of the elemental cube is 

pAxAyAz 

where p is the fluid density. Acceleration in the ^-direction may be written as the total 
derivative of the velocity 

Du . 

-X = UUx + VUy + WUs + Uf. 

Dt 

However, this value of the acceleration is only useful for a body of water in “absolute” 
coordinates. Our fluid exists in a rotating system on the surface of the earth, so a correction 
factor must be added to convert to a “relative” coordinate system. This correction factor is 
called the Coriolis acceleration and changes the net acceleration to 

UUx + VUy + WU~ + Uf — fv + fyW 

where / is the Coriolis parameter that is equal to two times the angular velocity of the earth 
times the sine of the latitude. The second Coriolis term is generally taken to be small for 
oceanographic models and we ignore it here. Equating the force terms with the product 
of the mass and acceleration terms yields the conservation of momentum equation in the 
rc-direction, which is 

p {UUx + VUy + VJU~ + Uf — fv) = — 


Px + (o'll)x + (&21 )y + (<73l)i 


( 3 ) 
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A similar derivation yields the y- and ^-momentum conservation laws. These equations have 
different terms for the Coriolis acceleration and the ^-momentum equation has a term for 
the hydrostatic effects of gravity. These equations are 

p (uvz + VVy + WV~ + Vf + fu ) = -Py + + (d-22 )y + (^32)* (4) 

p (uWx + vujy + wwf + Wi - fyU) = -pg -Pz + (<Ji 3 )i + (d 2 3)y + ( 033 )* (5) 

where g is the acceleration due to gravity. 

3.1 Governing Equations for the Two-Dimensional Model 

Equations (2)-(5) are the governing equations of our model. We have omitted all 
factors relating to temperature and salinity, two other properties that must be conserved in 
addition to mass and momentum. We also take density to be constant. For simplicity, we 
take p = 1 with no loss in generality; we do this by multiplying through the stresses and 
pressure by - p . 

In addition to the assumptions of the previous section, we will make another simpli¬ 
fying change. We first consider a set of governing equations with reduced complexity. We 
choose to ignore the ^/-direction and focus instead only on movement in two dimensions. Our 
goal is to analyze the impact of the horizontal coordinates relative to the depth coordinates. 
Assuming that all quantities are independent of the //-coordinate reduces our problem from 
a sloshing tank of fluid to a thin fluid sheet trapped between two panes on either side. In 
this context, the Coriolis acceleration has no meaning; the Coriolis terms come out of the 
governing equations when we remove all terms dependent on the //-coordinate. The four 
governing equations are reduced to three: mass continuity, conservation of momentum in 
the x direction, and conservation of momentum in the 5 direction; 

Ux + w z = 0 , ( 6 ) 

U£ + UUx + WUi + Px = (dll) £ + (dl 3 )i; * ( 7 ) 

Wi + UWx + WWz +pz = (<731)2 + (<5-33)2 ~ 9 , ( 8 ) 

where g is the acceleration of gravity. 

We assume the following constitutive relations among stresses and strains 

dn = 2 Eux, b -33 = 2Ewz , (7i3 = <5-31 = E{u & + u>x). (9) 

These relationships follow the form of the Navier-Stokes equations. However, the Eddy 
Viscosity (E) is substituted for the kinematic viscosity (u) because kinematic viscosity applies 
to laminar flow and is not appropriate for the turbulent flow of coastal ocean circulation. 
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+ X 

Figure 3: The domain of the two-dimensional model 


3.2 Boundary Conditions 

There are two surfaces which form the boundaries of the fluid in question. The lower 
surface is given by z = a(x) and the upper surface by z — a(x) + h(x,t ) where h(x,t ) is the 
height of the water column as seen in (Figure 3). We assume that fluid particles located on 
these surfaces at one time remain there for all times. The mathematical expression of this 
implication for the bottom surface is 

| (Z(i) - a(X :(i))) = 0 

using our Lagrangian notation. Noting the definitions in (1), we use the chain rule to 
differentiate the previous equation by t and express this boundary condition as 

w(x, a, t) — a,£u(x, a, t) — 0. (10) 

The complementary equation for the upper surface is 

o = | (z(f) - (a(A'(f)) + h(X(i),i))) 

= w(x, a + h,t) — axu(x, a + h,t) — hzu(x, a + h,t) — hf. (11) 

An additional boundary condition is that pressure at the free surface is the same when 
measured from either side of the boundary. We assume that the pressure of the atmosphere 
above the free surface is zero. Thus, the pressure at the free surface is also zero. This is 
expressed as 


p(x, a + h,t) = 0 . 


( 12 ) 
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Boundary conditions that equate the turbulent stresses on the upper and lower sur¬ 
faces to the frictional forces and wind forces on these surfaces are also required. We defer 
derivation of these conditions until we transform the physical equations into their non- 
dimensional forms because the introduction of scaling factors will greatly simplify that 
derivation. 

Lastly, we note that initial conditions are required for u, w , and h in order to solve 
this system of equations. 


Summary of Equations 


Mass Continuity 

Ux + w~ z = 0 

(6) 

^-Momentum 

iii + iiUx + wuz +px = (dn) £ + (di 3 ) £ 

m 

^-Momentum 

w t + uwx + wwz +Pz = (di 3 ) £ + (d 33 ) - - g 

(8) 

Lower Surface Boundary 

w(a ) — dxu(a) = 0 

(10) 

LIpper Surface Boundary 

w(x, a + h,t) — (a + h'j u(x, a + h, t) — h{ = 0 

(11) 
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4 Shallow Water Scaling 

We non-dimensionalize the variables in our equation by assigning a scaling factor to 
all dependent and independent variables. The independent scaling factors are 

x = Lx , 5 = Hz, t = jji, (13) 

while the dependent scaling factors are 

h{x,t) = Hh(x,t), a(x ) = Ha(x), 

u(x,z,t) = Uu(x,z,t), w(x, z, t) = ^-w(x, z, t), (14) 

p(x,z,t) = gHp(x, z,t). 

The scaling factor L is applied to horizontal coordinates while the scaling factor H is applied 
to vertical coordinates. U represents a typical horizontal velocity. These two scaling factors 
approximate the dimensions of the model’s domain. Since the hydrostatic pressure at a point 
depends on the weight of the water column above that point, the pressure is scaled by the 
product gH where g is the acceleration due to gravity. From these basic scaling factors we 
can see that x and u are related by the identity u = %. 


Following the convention of other works in this area [1], we define three dimensionless 
scalars A, e, and C 2 as follows 


A 


H , E r 2 _ gH 

L ’ t UL ’ U 2 ' 


(15) 


These scalars are called the aspect ratio (A), the inverse Reynolds number (e), and the 
inverse Froude number ( C 2 ). They appear as coefficients in our dimensionless equations. 
In most shallow water situations, the vertical dimension will be much smaller than the 
horizontal dimension. In the shallow water estuary problems, water depth is seldom over 50 
meters but horizontal distances are hundreds of kilometers. We therefore make the following 
assumptions about the relationships among these scaling factors: 


0 < A « 1, e«l, C 2 > 1. 


(16) 


The relationships between these scaling factors will be important for justifying mathematical 
simplifications of our governing equations. 


The scalings for the a terms follow directly from (13) and (14). These scalings are 


dn = 2 Eu £ 
d 33 = 2 Ew- Z 


EU_ 

2 r Ux , 


EU _ 
2^—w? 


L 


( 17 ) 
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£i 3 = <3- 3 i = E (u S + w&) 


EU f 1 


J u z + Aw s 


With scaling factors defined for all variables, the equations of motion and boundary condi¬ 
tions can be translated into dimensionless form. 


4.1 Dimensionless Mass Continuity Equation 

We non-dimensionalize the mass continuity equation (6) using the identities in (14) 
to obtain 


0 


Ux + Wz 
UUxXx + 


UH 

~~L 


W 2 Zz 


We evaluate Xj. and z z using (13) and reduce the above equation to 


or, 


U _ U _ 

—Ux + — Wz = 0, 


U S +Wz = 0 . 


(18) 


4.2 Dimensionless ^-Momentum Equation and the Hydrostatic 
Approximation 

Using the scaling factors from (13)-(17), the z-momentum equation (8) becomes 

U 2 H U 2 H . . U 2 H , 2 . EU (l \ EU , . 

+ -jy (uw ^ + ~U~ ^ ^ + 9P ' Z = ^ X u ~ z + Xw *) _ + ( w *)z - 9 

which can be written as 

A 2 [-iDf + (uw)- + (w 2 ) _] + C 2 (pz + 1) = e [(% + X 2 Wx) _ + 2 (wg)J[ (19) 

Until this point our model has been exact; no terms have been dropped to simplify 
the problem. Now with the z-momentum equation (19) and the scalars in (15), we have 
the opportunity to do so. In the z-momentum equation all terms are multiplied by one of 
the factors A, e, or C 2 . Recalling the relationships between scaling factors from (16)—that 
0 < A < 1, 0 < e < 1, and C 2 3> 1—we see that the other terms of the z-momentum 
equation are insignificant compared to the term multiplied by C 2 . Removing these, our 
approximation is 

C 2 (p- Z + 1) = 0. 

Now we integrate this last identity with respect to the vertical coordinate from an arbitrary 
point 7 to the free surface of the fluid at a + h to get 

/ a-\-h 

(pz + 1) dp = p(x , a + h,t) — p(x , z,t) + (a + h) — z = 0. 
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Since the pressure at the free surface is zero, we obtain 

p = a + h — z. (20) 

This last identity is referred to as the hydrostatic approximation. 


4.3 Dimensionless x-Momentum Equation 

The ^-momentum equation non-dimensionalizes similarly. From (7) we have 

A + {u 2 )~ + (wu) s +px = (<7h) 4 + {a 13 )~ 

We replace the scaled parameters with the non-dimensional versions to get 


U 2 _ u \_ 2 v u 2 , gH _ ( 

T Ul + T + T ^ wu ^ + ~l p * = 2 ( 


EU 

U 


u x + 


EU (l 

HL \ UI + XW ‘ 


or 


Uf + (u 2 ) _ + (wu)z + C 2 p s = 2 (eu.x)^ + _ + ( ew x)- z 


with C 2 , e, and A as defined in (16). 

In the situation where e is constant, we exploit the identity w? = — % (18) to reduce 
the right hand side of the last equation. Also, we use the hydrostatic approximation (20) to 
replace p s with a s + h s to obtain 


ui + (u 2 ) _ + {wu)- + C 2 (% + h s ) = (eu, s ) s + (j 2 U s) _ 


( 21 ) 


4.4 

to 


and 


Dimensionless Boundary Conditions 

The boundary condition at the bottom surface (10) and the free surface (11) transform 


w (. x , a(x)) = 

(22) 

(x, a(x) + h(x,t),t) = hf + u(a + h)x- 

(23) 


4.5 Dimensionless Surface Traction Approximation 

The derivation of the surface traction approximations in (A) is taken from [4], The results 
of this derivation are two boundary conditions, (A-8) and (A-9). These can be written as 


^^x^x 


= —K {r hu(x,a,i) 

A z 


(24) 
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for the lower surface z = a and 

- e (as + h s ) u s = y K ' v ( Mw _ 5 + t)) (25) 

for the upper surface where z = a + h. K w is the coefficient of friction between the fluid and 
the air at the free surface, Kf r is the coefficient of friction between the fluid and the bottom 
surface, and u w (x,t ) is the velocity of the wind at the free surface. 


Summary of Equations 


Mass Continuity 

Ux + Wz = 0 


(18) 

x-Momentum 

% + (« 2 )s + + C 2 (a s + hx) 

= (ehs)s+ (^h 2 -)_ 

(21) 

Upper Surface Boundary 

w (x, a,t) = ueix 


(22) 

Lower Surface Boundary 

w (x,a + h, f) — ht + u(a + h)x 


(23) 

Lower Boundary Traction 

jiUz ~ eax.Ux = jiK b hu(x , a, t) 


(24) 

Upper Boundary Traction 

~\2^z ^ (aT h-s) its ^2 L\ v (u w 

— u(x, a + h, f)) 

(25) 
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5 s -Coordinate Transformation 


We now introduce this fundamental transformation to map our equations of motion 
from a variable domain in the vertical direction to a fixed domain in the vertical direction. 
This simplifies our bathymetry and moves complexity from the domain to the governing 
equations by replacing a formerly irregular domain with a regular shape with constant height 
1. We let 


z — a(x) . . . . T ._ 

x = x, t — t, s = T f _ , a{x) = a{x), h{t,x) — h{x,t). 


Then 


and 


h(x, t) 


rv r> &X T Shx r> rv 1 n rv rv Silt r> 

dx — d x - - O s , Oz = yd s , Ot — O t 7—O s , 

n n n 


(26) 


(27) 


u(x,s,t) — u(x, z,t), w(x,s,t) — w(x, z,t). (28) 

For the free boundary of the fluid domain—the interface between wet and dry areas where 
h = 0—division by h produces a singularity. This problem will be addressed in later chapters 
as we numerically approximate the governing equations. Note that on the bottom surface 
z = a(x) and 

a(x) — a(x) 

^5 


s = 


h(t, x) 


whereas on the free surface z = a{x) + h(x,t ) and 

(a(x) + h(x,t )) — a(x) 


s = 


h(x, t) 


= 1 . 


We have changed the physical domain from a region defined by a(x) < z < a(x) + h(x,t) 
to a vertical domain defined by 0 < s < 1. This simplification comes at the cost of adding 
complexity to the governing equations. 


Following the precedent of previous works [1], we introduce the function 0, defined 
by 

4>(x, s,t) = w — (a x + sh x )u — sh t 
or 

w = (f) + (a x + sh x )u + sh t . (29) 

The function cj) has no simple physical interpretation. Instead, it is a mathematical artifact 
that exists to simplify the governing equations by removing w. In both the mass continuity 
and re-momentum equations, all terms containing w can be replaced with an equivalent term 
containing <fi. An important property for simplification is that (j) satisfies the boundary 
conditions 

4>(x, o, t) = 4>(x, l, t) — o 

which can be seen by putting (22) and (23) into (29). 


( 30 ) 
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5.1 Mass Continuity Equation in s-Coordinates 

The mass continuity equation (18) transforms into s-coordinates as follows: 

0 = Ux + Wz 


( cl x T sh x 

= u ’-{—i^) u ’ + h w - 

If we multiply the last equation through by h we obtain 

0 hu x (ci x T sh x )u s T uq, 

Differentiating the function 0 defined in (29) with respect to s yields 

w s = 4>s + h x u + (a x + sh x )u s + h t . 

If we replace w s in (31) by (32), we see that the continuity equation reduces to 

ht + ( hu) x + 4> s — 0. 

This equation is exact. 


(31) 


(32) 


(33) 


5.2 ^-Momentum Equation in s-Coordinates 

When the ^-momentum equation (21) is expressed in s-coordinates and multiplied 
through by h, the left side of this equation becomes 

( hu t — sh t u s ) + ( huu x — (a x + sh x )uu s + huu x — (a x + sh x )uu s ) 

+ (wu s + w s u ) + C 2 h (a x + h x ) = ... 

Noting that a s and h s are both zero allows us to remove those terms. We now use w and w s 
as defined in (29) and (32) to simplify the left hand side to 

(hu) t + ( lm 2 ) x + (0u) s + C 2 li (a x + h x ) = ... 

Applying the definitions in (26)-(28), the right hand side of the r-niomentum equation 
(21) transforms to 


••• = (ew*) s + T o U z 


e 

A 2 
a T + sh 


I ^ i ^x \ ^X + sh x ( OLx shx \ ( £ 1 ( 1 

= 1 5— “-j. - —5— ( £ “* “ £ —— + (vS [h u ‘ 

If we multiply this expression through by h and note that we can move h in or out of an s 
derivative because h is independent of s, we find that the right hand side reduces to 


e \ , / / a x + sh x 

u s ) +h[e(u x - 7 - u s 


X 2 h 


h 


. a x T sh x 
£ \^x sh x ) ( 'U'x ^s 
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Summary of Equations 
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6 Separating Mean and Residual Velocities 

The variable u(x, s, t) represents the lateral velocity of the fluid in the artihcal domain 
created by the s-transformation. The x-direction is the only direction of free motion in this 
domain, while in the physical domain there is motion in both the x- and ^-directions. The 
horizontal velocity varies over the water column of the fluid; this is expressed mathematically 
through the dependence of u on s. We split u into two components: one that is depth- 
dependent, and one that is independent of s. This decomposition is expressed by 

u(x,s,t) — u(x,t) + u(x,s,t) (37) 

where 

u(x,t) — / u(x,s,t)ds, 

Jo 

u(x , s, t ) = u(x, s, t ) — u(x, t). 

We refer to u as the depth-averaged horizontal velocity and u is the horizontal velocity 
residual. Additionally, u satisfies 

u(x, s , t) ds = 0. (38) 

6.1 Depth-Averaged Mass Continuity Equation 

When we apply the decomposition (37) for u to the mass continuity equation (33), 
the result is 

h t + ( h(u + u)) x + (f) s = 0 . 

If we integrate this equation over s from 0 to 1 and exploit the boundary conditions in (30) 
and the identity (38), we obtain 

ht + (hu) x = 0 (39) 

With this transformation, the mathematical artifact <fi is removed from the mass continuity 
equation. This was the purpose of replacing terms dependent on w with (f>. While calculation 
of values of w is not important to determining horizontal velocity fields, we can still use 
information about (p to determine values for w in post-processing. If we subtract (39) from 
(33) we see that 

(ft s (hu ) x • 

We then use (f>(x, 0 ,t) = 0, one of the boundary conditions of (p in (30), to obtain 

(p(x, s,t) = — J u(x,r},t) d'q'j . (40) 

Because of (38), the function <p given in (40) satisfies <p(x,l,t) = 0, the other boundary 
condition in (30). We now have an expression that we can evaluate to solve for cp and 
by extension w. Although results for these values are not calculated in this project, that 
calculation would be a useful extension of this work. 
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6.2 Depth-Averaged x-Momentum Equation 

We apply the decomposition in (37) to the ^-momentum equation (34) to obtain 

(h(u +u)) t + (h(h 2 + 2 uu + u 2 )) x + (0(h + u)) s + C 2 h(a x + h x ) = + “)«) + 


he ( (u + w) x - 0,1 (u + h) 



e (a x + sh x ) (w + u) x - Q ' x (u + «). 


Separating terms and noting that u s = 0 yields 


(. hu) t + (hw) t + ( hu 2 ) x + 2 (huu) x + ( hu 2 ) x + 0 s w + (0h) s + C 2 h(a x + h x ) = 


A 2 h 


u s ) + 


t / /= “I - shx ~ 

he ( (-u + u) x 1 —-—-w. 



/ 7 \ / /= “1“ s/i# - 

e (a x + sh x ) ( (u + u) x ---w. 


We now integrate the last equation with respect to s from 0 to 1. Noting that 
Jq 1 uds = 0, /q 1 hds = h, J^ uds = u, and that 0 is zero at both s = 0 and s — 1, we find 
that the left-hand side of the integrated expression is 


({hu) t + (hu 2 ) x + 



+ C 2 h(a x + h x ) — ... 


Integrating the right hand side we obtain 


e 

A 2 


e 

A 2 


h s (l) u s ( 0) 


h 


h 


+ ( heu x ) - 


e (a x + sh x ) ds 


/ , , \ I /= . flx + sh x _ 

e (a x + sh x ) ( (w + w) x - - -u. 


u s (l) « s (0) 


+ (heu a 


— e 


e (a x + h x ) -u(l) - ea x w(0) 
cl x “h h. 


h h 

_i i- 

(a x + h x ) ((u + u) x (1) - 0,1 + hx u s ( 1) j - a x ((ft + u) x (0) - ~u ( 


h 


J^u s ( 1) - e(a x + h x ) (w x (l) - a ' x + /7a ’ u s (l)^ - ^u s (0) + ea x (u x (0) - 


+ {heu x 


^9 7i w ( , u w u "u(l)) hK{ T (u -\- w(0)) 4- (heu x ) x 


e (a x + h x ) «(1) - ea x h(0) 
e (a x + h x ) h(l) - ea x h(0) 
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Thus, the depth-averaged ^-momentum equation is 


(hu)t + {hu )# hu ds J C h{u x hx) — 

-^K w (u w -u- 71(1)) - ^hK ir (u + 7t(0)) + ( heu x ) 3 


6.3 Residual x-Momentum Equation 


e (a x + h x ) u( 1) - ea x u( 0) 


(41) 


To obtain the equation for u, we subtract the depth-averaged ^-momentum equation 
(41) from (34). The result is 

(■ hu) t = - 2 (huu) x - ( hu 2 ) x + (^ j hu 2 ds) - (j) s h - (<jm) s - ^ [AT* (u w - u - 77(1))] 


[~hKg (■u + 77(0))] + [e (a x + h x ) 77(1) - ea x u( 0)] + ^ (u s ) s + ( heu x ) x - 

// 7 \~\ /^/ 7 \ = ~\ &X Shx ~ 

(e (a x + sh x ) u s ) x - I e (a* + sh x ) I [u + u) x --- u t 


(42) 


We lot 

define the function G by 

G = — h t u — 2 (. huu) x — 


D A 2 ’ 


h (u 2 ) — h I u 2 (x, r /, f) dr] 


- (p s u - (4m)z 


+ [e (a x + h x ) u( 1) - ea x 7i(0)] + (heu x ) x - (e (a x + s/i x ) 7t s ) x . 
and note that we have selected the terms for G from (42) so that 

/ G(x, s, i) ds = 0. 

Jo 

The equation for u then becomes 

hu t G T —j- [^1 T A (oj X sh x ) ^ [c (u x T sh x ) (u x T n x )]^ 

— D [A' w (n w — u — 77(1)) — hK {r (u + 77(0))]. 

6.4 Equations for Boundary Conditions 

On the lower surface s — 0, (35) becomes 

Y (l + A 2 a 2 x ) 77 s (0) — ( ea x (u + 77(0))) = DhK(u + 77(0)) 


(43) 


(44) 



u-u( 1)) (45) 


and on the upper surface at s — 1, (36) becomes 

— (l + A 2 (a x + h x ) ) « a (l) — (e (a x + h x ) (u + ft(l))) = DK ® (n w — 
In the following section, we make the assumption that 


6.5 Approximating the x -Velocity Residual 

The derivation of the approximate equation for u follows the approach of Greenberg 
in [3]. The approximate equation is obtained by solving the velocity residual (43) with the 
boundary conditions (44) and (45) and exploiting the fact that D is much larger than both 
e and 1 and that D A 2 is much smaller than D. We retain only terms on the right hand side 

of (43), (44), and (45) that are 0(D) while ignoring terms that are 0(1), 0(e), and 0(D\ 2 ). 

Thus, we investigate the following system of equations 

hu t = ^u ss - D [K w (u w - u - h(l)) - hK {r (u + u( 0))], (46) 

w(0) = h 2 K fr (■ u + ft(0)), (47) 

u s (l) = hK w (u w - U - u(l)). (48) 

The solution to (46)-(48) may be written as 

u—U+-w 


where U is an approximate steady state solution satisfying 

U ss = (kJi(u w - u) - K tr h 2 u - K w hU( 1) - K ir h 2 U(0)^j , (49) 

U a {0) = K fr h 2 (S+L7(0)) t (50) 

U s (l) = K w h(u w ~n-U(l)^j (51) 


and w is the transient component which may be written as 

OO 

w(x,t) = WkeO fJ ’ kt ' > (p k (s ) 

i =1 

where <p k and /i^ satisfy the eigenvalue problem 

p (W), («) + A»fc/(1) + K t hV( 0)) + n^(s) = 0, 
0‘(O) = K<,h 2 4>\ 0), 

^( 1 ) = 1 ). 


0 < s < 1 
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These eigen-functions satisfy 

(p k (s)ds = 0, k — 1,2,3,... 

0 fc ( s )0" , (' s ) ds = 0. k — 1,2,3,... 

If we normalize 0 fc (s) to satisfy 

(<p k ) 2 ( s) ds = 1 

we find that 

0 < Wt = ^ /’ Of W * + § (x w h (/(l)) 2 + A'f,/! 2 (/(o)) 2 ) . 

Since the only solution to 

Ms + (AT w /i0(l) + /i fr fi 2 0(O)) = 0, 

</> s ( 0) = K fl h 2 (/>( 0), 

&(1) = -Aw^(l) 

is 0(s) = 0, we are guaranteed that n = 0 is not an eigen-value. Our assumption that 
D = e/X 2 S> 1 guarantees that solutions to the transient problem decay rapidly for any 
initial condition. This allows us to use the steady state solution U as an approximation 
for u. Therefore, terms in the depth-averaged ^-momentum equation that contain u can be 
replaced with U. 




The solution to the steady-state component of the residual (49)-(51) may be written 


as 

2 

U = U( 0) + A fr + f/(0)) s + (K v h («w - u) - K ir h 2 U - K w hU( 1) - K h h 2 U( 0)) y (52) 


where 


U( 0 ) 

m 


B 


-K w h(u w -u) (1 + K w h/4) K b h 2 u 
6B 6B ’ 

(1/3 + Af r fi 2 /12) K w h (u w — u) Kf r h 2 u 
B + 6 B 

, (K w h + K {r h 2 ) K w K {r h 3 

' n ' -t i-v 


(53) 

(54) 

(55) 


and 
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6.6 Summary of Depth-Averaged Equations 

There are two depth-averaged equations of motion in our problem, summarized in 
the box below. Recall that U is the steady state approximation to u which we will use in 
subsequent computations. U does not depend on x or t directly, instead it depends on u, 
u w , and h, and those values depend on x and t. In the next section we solve the two partial 
differential equations of motion using numerical methods. 
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7 Finite-Dimensional Approximating Equations 

Until now, we have developed a continuum description of the fluid flow where all 
fields, particularly u and h , are defined at all points x and at all times t where the water 
column h(x, t ) is positive. This continuum problem is infinite dimensional—for example, a 
Fourier series representation of the solution will, in general, need infinitely many Fourier 
coefficients. Because the system formed by (39), (41), and (52)-(55) is non-linear, obtaining 
an exact solution is impossible. Therefore, this section presents a method for obtaining an 
approximate solution based on a finite-dimensional approximation of the infinite-dimensional 
problem. 

We use a Lagrangian approach to create this approximate solution. With this method 
we define a series of “fluid-marker points” that we track for all time. In the context of this 
model, the term fluid-marker point does not refer to any specific particle at any place in the 
fluid. Instead it refers to the projection of a column of fluid onto the model’s “grid”. The 
grid is a structure for describing the lateral space of the fluid domain. In the two-dimensional 
model, the lateral space of the domain has one degree of freedom and the fluid-marker points 
are projections of a vertical column of water onto this line of lateral space. 


We assume the wetted fluid region of the continuum problem has a lower bound of a 
and an upper bound of b at t — 0 so that h(x, 0) satisfies the following constraints 

0, x < a, 

h(x , 0) = h°{x) >0, a < x < b, 

0, b < x. 


We create a grid that divides the region between a and b into N + 1 points. If this grid is 
uniformly spaced then the i th point is denoted by 

x — 1 

x° = a H-— — (6 — a), 1 < i < iV + 1. 

i N 

\: <: all grids will be uniformly spaced; depending on the initial conditions, we may find 
it preferable to distribute grid points by other means. However, no matter the method of 
dividing the lateral space, there will always be N + 1 grid points that are identified as fluid- 
marker points. The endpoints of this grid match the edges of the continuum problem where 
x° = a and x° N+1 = b. The difference between two grid points at t — 0 is 

- A = 1 <i<N. 

Given the sequence of grid points and the function for the initial water column 
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h(x, 0), we define the sequence {Mf}^ 1 by 


Mf 


0, i = 1, 

f*J h°(x) dx, 2 < i < N + 1. 


The numbers M° represent the amount of fluid in the interval [x?, xf ] at time t — 0 which 
we refer to as “fluid elements”. We compute the amount of fluid in each interval [x°,x° +1 ) 
by 

< = M? +1 - Mf, 1 < i < iV 

and the average height of the water in the same interval by 


K = 


mx 


N 


-rn. 


x] 


i +1 


X; 


(56) 


Now that the height of the water column in each interval is approximated by the average 
water column, we no longer have any singularities in our governing equations along the free 
boundary where h = 0. 


Given any domain of fluid described by a continuum mechanics model at a time t = 0, 
we have a finite-dimensional or discrete approximation of that domain. Our discretization 
will be a set of ordinary differential equations that tracks through time the position of the 
grid points the velocities of these fluid-marker points, and the average height of the 

water column between successive marker points. We adopt the following notation to denote 
these three values: 

(1) Xi(t) is the position at time t of the marker point that was located at x° at t — 0 
for 1 < i < N + 1. 

(2) Ui(t) is the velocity of the marker point X;(t) at time t for 1 < % < iV + 1 and 
represents the value of u(xi(t),t). 

(3) hi(t) is the average height of the water column in the interval [&*(£), Xi + i(t)) for 
1 <i< N. 


Note that x and u satisfy 


C LJu 1 - , , , , 

— = Ui(t), 1 < I < N, 

Xi( 0) = 1 < i < N. 


and 
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7.1 Ordinary Differential Equation Approximation 

Our strategy is to replace the partial differential equations in (39) and (41) by an 
approximating system of 2N + 2 ordinary differential equations. At each of the points x t , 
from i = 1 to % = N + 1, we create a pair of ordinary differential equations of the form 


d,Xi 

dt 


= u h 


and 


<-1 


d‘ u i~ i „ / o , o 

+ 3 (' m i-i + m 


\ d u i 0 du i+ i 


F(xi,Ui), 


where m; is a constant mass coefficient and F(xi,Ui ) represents forcing terms. This creates 
N + 1 pairs of ordinary differential equations, each corresponding to the different index 
values of i. All 2N + 2 equations are expressed more succinctly in vector form. We define 
two vectors 

x = {^}f= + i 1 , (57) 

and 

u = {«,}£' (58) 

and rewrite our approximating ordinary differential equations as 


dx 

~dt =U ’ 

- du _ 

M -= F( X, U ) 


(59) 

(60) 


where M is a tri-diagonal matrix and F is a vector of forcing terms. We calculate M and F 
in the following sections. 


7.2 Discretizing the Depth-Averaged Mass Continuity Equation 

We integrate the continuum form of the averaged mass continuity equation (39) over 
the interval [xi,Xi + 1 ]. Assuming that all instances of h t are evaluated at t, the result is 

rx i+ 1 (t) rx i+ i(t) 

0 = / h t dx+ (hu) x dx 

J Xi(t) J Xi(t) 

rf r x i+i(t) 

= — / hdx - (x i+ i) t hi(x i+ i) + (xi) t hi(xi) + u(x i+1 , t)hi(x i+1 ) - u(x i ,t)h i (x i ). 

at Jxi(t) 


Since (xi(t)) t is equal to u(xi(t),t ) the above expression reduces to 


d r x i+i(t) 

— hdx = 0. 

dt Jxi(t) 


(61) 



34 


The integral of the water column over an interval between fluid markers is the amount of 
fluid in that interval, or, since we assume that density is a constant p — 1, the mass of that 
interval . From (61) we see that this mass does not change over time so we generalize that 
m° = m,. We see that 



h(x,t ) dx 



h(x, 0) dx 


f x i +1 

/ h°(x ) dx = m° 

Jx? 


and (56) implies that 


mi = (x i+ i(t) - Xi(t )) hi(t). 


This is the discrete approximation to the depth-averaged continuity equation. 


(62) 


7.3 Discretizing the Depth-Averaged x-Momentum Equation 

We have established a grid scheme that discretizes the fluid domain. Now we derive 
a set of evolution equations for the sequence of velocities {-Uj(f)}^ 1 which are consistent 
with the depth-averaged x-momentum equation. To obtain these equations we extend the 
discrete velocities (Sj(f)}^ 1 and average water columns {hi}f =1 to fields defined on the 
interval (—oo, oo). For iq, we do this by piecewise linear extension: 

{ Ui(t), x < Xi(t), 

Ui(t) + i x - x i(t))i x i(t) < x < x i+ i(t) and 1 < i < N, (63) 

u N+ i(t), x N +i(t) < x. 

For the average water columns, we do this by a piecewise constant function 

0, x < xi(t), 

hi(t), Xi(t ) < x < x i+ i(t) and 1 < i < N, (64) 

0, x N+ i(t) < x. 


h(x, t ) = 


where hi(t) = 


Xi+l(t)-Xi(t) ' 


To obtain the evolution equations from the discrete velocities {■Uj(f)}^ 1 , we integrate 
the depth-averaged x-momentum equation around each point x*. We choose to integrate from 
the midpoints between adjacent fluid markers. We let 


Xi- 1/2 


Xj+Xj-t 

2 


and x i+ i/ 2 


2 ’ 


be the lower and upper bounds of this region of integration and a corresponding notation 
for u is 
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The depth-averaged x-rriornentum equation (41) can be reorganized as 


momentum 


pressure 


c 2 

( hu) t + ( hU 2 ) x = - hU 2 ds ) - — (h 2 ) x - C 2 ha 

+ ~^K w (u w -u-U{ 1)) - ^hK {r (u + U(0)) + (heu x ) x - 

S. 


potential 

^2? 


e (a x + h x ) U (1) - ea x U (0) 


wind and friction 


viscosity 


remainder 


For simplicity, we choose to ignore the remainder term because the remainder term is much 
smaller than the other terms in the depth-averaged x-niomentum equation. We integrate 
around each point ay to obtain 


momentum 


c i+1/2 


(( hu) t + ( hu 2 ) x ) dx = — C 2 


' x i- 1/2 




+ 


(ehu x ) x dx + 


v i~ 1/2 


potential 

r x i+l /2 ' 

/ lia x dx — 

d x i-l/2 
f X i+ 1/2 e 

- , A2 

v i- 1/2 


pressure 


*+1/2 


i—1/2 L 


CM ? ) + h I U' z ds 


dx 


K w (u w -U- U{ 1)) - hK b (u + *7(0)) 


da; 


viscosity 


wind and friction 


(65) 


7.3.1 Momentum Terms 

The momentum terms from the integrated averaged x-niomentum equation (65) can 
be written as 

(hu) dx - hiU i+ 1/2 ( 24 + 1 / 2 ) t + (xj_i /2 ) t + hi (u i+1 / 2 ) 2 - d*_i (i*-i/ 2 ) 2 • 


d 

dt 


r x i + 1/2 

E i-1/2 


Recalling that ^ — u, this last expression reduces to 


d f Xi +^ 

dt J, r , , o 


(hu) dx. 


( 66 ) 


Using (62)-(64) we see that 


K i + l/2 


1/2 


(dw) da; = hi 

d/ 

1 


c i + l/2 


u dx + hi-\ 

Uj + l+Uj 


u dx 


Jx i ~ 1/2 

2 - \ ay (1 - a;* 


+ hi -1 


Wi + 


Uj+Uj-i 




+ 3aq (mj_i + mi) + aq+im* 
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Therefore, (66) can be rewritten in vector form as 


, „ . dn 

M 'n 


where M is a tri-diagonal mass matrix with elements defined by 


M L i = -mi, 

O 

3 

Mjv+i,jv+i = -m n, 

O 

3 

m j,3 = i + '"./)• 

m jj-i = l ni .i i- 

M ij+i = 


2<J<N, 
2<j<N + l, 
1 <j<N- 


(67) 


7.3.2 Potential and Pressure Terms 

The potential term from the integrated, depth-averaged x-niomentum equation (65) 
can be written as 

C 2 hi -1 (a(xi) - a(xi- 1 / 2 )) + C 2 hi (a(x i+ 1 / 2 ) - a(ay)) . (68) 

The pressure terms are 

%~h 2 i ~ %-h 2 ^ + hi [ U{x i+ i/ 2 ) 2 ds — /q_i [ U(xi_ 1/2 ) 2 ds. (69) 

* * Jo J o 

These terms can be evaluated in the forms above. We represent their vector forms by F pre s 
and F p Q t. 

7.3.3 Viscosity Term 

In evaluating the viscosity term, we assume that e is piecewise constant in the intervals 
between fluid markers, so that 

e(x, t ) = 6j, Xi < x < x i+ i and 1 < i < N. 

Using this e, we write the integrated viscosity term as 

7 ^h+1 ^ l 7 Hi Ui —l 

Xi+1 %i— 1 

There are many possible choices for e t . One choice is 

2 


e,; = 


H ( b — a 


8 dt V N 


(70) 
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This value for e has the advantage of being constant in both space and time. Note that dtt 
is the size of the time step that we will be taking in our numerical computations. Another 
choice is 




(x i+1 - XiY 


8 dt 


The viscous forces may be written in vector form as 

Mu 

where u is defined in (58) and M is a tri-diagonal matrix defined by 

e x hi 


(71) 


(72) 


M,, = 


Aaq ’ 


M 


N+1,N -|-1 


Ax 


N 


M 
M 

We are using the notation 



rb 

1 

1 

3,3 ~ 

V Ax 3~ 1 



3,3- 1 

< 

1 



iJ+i 

AxY 


+ 


AXj 


2<j<N, 

1 < j < N. 


(73) 


Axi = (x i+ 1 - Xi ). 


7.3.4 Friction and Wind Terms 

The wind term from the integrated averaged x-rnornenturn equation (65) is 

r x i+ 1/2 e 


A 2 


K w ( u w — u — 7/(1)) dx 

Jxi-yi L V y J 

and the friction term is 

r x i+i/i c r / _ \ ' 

—hK b (u + f/(0)J dx. 

Using the definitions of t/(0) and f/( 1) from (53)-(55), we expand these two terms to 

dx 


x i+l/2 g 
7-7/2 ^ 


i+1/2 ^ 

A 2 


and 


x i — 1/2 

f x i+ 1/2 g 

~\2 

Hi-l/2 


r , , =, 4 + K fv Klh\. K fr K w h ' 

A w («w - u) - 7777 ; h (u w - U) -tv; hu 


12 B 


6/5 


„ ,= , K w K{ r h . (1 + K w h/ 4) Kf h 2 = 

-K{ r hu H ——h (u w - u) H-tv;- hu 


6B 


6B 


dx. 
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Using (62), we evaluate these expressions to obtain the combined wind and friction term 
K w 6i-i f 4 + (K w hi -1 + 2) K{ r hi-i 


(«W - u)i_ 1 


8A 2 


hi -1 


, -n 3 AT, 

+ (u w - u) 


12i3 

4 + ( K w hi-\ + 2) Af r hj_i . 

U—i 1 hi—i I Axf— i 


e,: 1 


8A 2 

4 + ( K w hi + 2) A fr A 


125 


125 


hi Axi 


+ («w - «)i 


i+1 


K e 

1 v W c 2 

"8A 2 ” 


4 + (K w hi + 2) K{ r hi . 
125 


+ 


Ad r e*-i / (1 + K w hi-i/A) Kfrhi-! - K,, 


Ui 


L 8A 2 V 

3Af r 


65 


hj_i — 1 m,_ 


8A 2 


U-i 


/(I + AAvh,;_i/4) AfA_i — A w 

v 65 


hi— i - 1 ) 


+U 


/ (1 + K w hj/ 4) AfA — 


■hi — 1 ) tri 


+ Ui+l 


65 

Kfr€i / (1 + K w hi/ 4) Kfrhi - AU 

[ 8A 2 V 65 


hi — 11m. 


To express this expression in vector form, we define two matrices Mf and M w . The matrices 
are 


(M 


w) 1,1 


K w ei ( 4 + (K w hi + 2) Kf r h\ 


8A 2 


1 - 


125 


hi Ax 1} 


(M w )jv+i,jv+i — 


K w £n f i 4 + ( K w h]y + 2) Kf r h 
8A 2 


1 - 


v 


(M 




3 A\ v 


U-i ( 1 — 


125 

4 + (K w hi-i + 2) Khhi-i 


(M v 


“ 8A 2 

(M w )jj + i 


125 

4 + (K w hi + 2) AfiA 
125 

Ah-i A 4 + (K w hi-i + 2) Kf r hi-\ 


hNJ Axjy, 

hi-1) Axi-f 


+Q I 1 


K ) Aon 


2 < j < N, 


I< e 

1 v W C Z 

"8A 2 ” 


125 

4 + (Kwhi + 2) Kfrhi 
125 


hj_i Ax*-!, 2 < j < A + 1, 


hi A Xi, 


1 < j < A, 
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and 


(M 


frJU gA 2 

(M 


K h e i /(1 + K w hi/ 4) K tr h, - AT, 


6 B 


hi - 1 m i, 


fry N-\-l,N-\-l 

3Af r 


KfrtN /(I + K w Hn/A) K^Hn — A" w 


6£> 


/ijv - 1 mjv, 


8A 2 


Q-i 


(1 + A^ w /ij_i/4) K{ r hi-i — A w 


6£ 


/i,-i - 1 mt-i 


+Q 


(1 + K w hi/A ) K ir hi - K\ 

as 


hi - 1 m, 


(M 

(M 


frA'd-i — 


Af r ej_i / (1 + K w hi-i/A) Kfrhi-i — K v 


8A 2 


6£ 


2<j<N, 


K-1 - 1 mi-i, 2 < j < N + 1, 


fr/jj+i 


8A 2 

We also define a vector 


Af r e,; / (1 + K w hi/A) K t - r h t — A w 

as 


hi - 1 mi, 


1 < j < N. 


u w = {(MwWj 1 • 


With these definitions, the friction and wind terms may be written in vector form as 

M w (u w - u) - M fr u. (74) 


7.3.5 Summary of All Terms of Depth-Averaged s-Momentum Equation 

When the terms of (65) have been translated into vector form, the equation becomes 

pressure potential viscosity wind and friction 

Fpres - Fpot + Mu + M w (u w - u) - M fl u . (75) 

This is the equation in the form of (60) that forms the second part of our ordinary differential 
equation approximation, the first part of which is 

dx 

— = u. 
dt 

We solve this system of ordinary differential equations through operator splitting. Each time 
we step forward through time by one time step, we evaluate our system in two distinct parts. 
For the first part we hold x constant and update u; in the second part we hold u constant and 
update x. We use a time grid where the constant interval dt is the length of time between 
t n and t n+1 for any n. Another way of stating this is that t n = ndt for n = 0,1,... etc. 
We use a superscript notation to signify the time step at which a value is being calculated, 


momentum 
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thus u n+1 refers to the depth-averaged velocity one time step after u". 
be described by 


Step 1 

Step 2 

t n <t< t n+1 

t n <t< t n+1 

f = o 

— = u 

dt 

M|=f(x,u) 

= 0 

at 


The two steps can 


(76) 


Note that we do not consider u w a variable upon which (75) depends; it is given or externally 
supplied data like the mass, viscosity, and friction matrices. 
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8 Implementation of the Two-Dimensional Model 

At the end of the previous section, we found a solution for the velocity fields of the 
two-dimensional model. This solution involves a “leap-frog” updating process and uses a 
Lagrangian method to track particles from time zero through all future time steps. Within 
each time step, we solve first for the depth averaged mean velocity, u, of the fluid at each 
particle location using data from the previous time step. Next we use the mean velocity to 
update the location of each particle, x. Finally, with a new array of particle locations and 
a constant array of masses, we can calculate the water column, h, in each interval between 
particles. 

There are many options for implementing the solution of the final differential equation 
in MATLAB. Two of the easiest solutions are the explicit and implicit Euler methods. In 
explicit methods, such as the forward Euler method, the value of the variable of interest—in 
this case the depth averaged velocity, u —at the next time step is defined in terms of known 
values from the current time step. Implicit methods, such as the backwards Euler method, 
defines the variable at the next time step in terms of unknown values from the next time step. 
The variable of interest must then be solved for; this solution takes time but the advantage 
is that implicit solutions are usually more stable allowing the model to use larger time steps. 
The two-dimensional model has been implemented using both approaches. 

The final differential equation (75) represents a series of n equations, one for each 
tracked particle. In the context of implementing this model, the advantage of the explicit 
Euler method is that with certain viscosity values each equation can be independent. That 
means that there is no need for linear algebra when solving; instead MATLAB will just 
solve a series of n independent equations. The implicit Euler method does require linear 
algebra because all the differential equations in the series are dependent on each other; but 
the advantage of this method is its increased stability. The implicit method produces stable 
solutions for almost all initial conditions when the ratio A is 0.1 while the explicit method 
can fail if the spacing between the particle locations, dx, becomes too small. Testing of the 
two methods against each other have shown that the explicit and implicit methods provide 
near identical results when both are stable and initial conditions are the same. However, 
when the residual velocities and wind and surface friction forces are added, the domain of 
inputs on which the explicit solution is stable becomes much smaller unless the ratio ^ is 
reduced, so the implicit method is preferred. 

Whether the explicit or implicit method is used, MATLAB implementation of this 
leapfrog process takes as input the state of a region of water at any time and produces 
as output the state of that water at a future time. Here state is defined primarily by the 
three variables of interest to this program: depth averaged mean velocity at each particle, 
u , location of each particle, x, and water column of each fluid element, h. In order to 
perform this calculation, the program needs these three pieces of information at time zero. 
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In addition, the program needs a variety of other initial data, including a function defining 
bathymetry, a function defining wind, constants for wind friction and bottom friction, initial 
velocity of the fluid, and more. 

The program provides two graphical outputs for visualizing results of the model. A 
plot of height over bathymetry versus x gives the user a view of what the two dimensional 
model might look like in a real physical setting-or at least as real as one can get with two 
dimensions. The other graph plots the center of mass versus the average momentum of the 
mass of water. This graphic is useful as a diagnostic. For cases where the bathymetry is a 
parabola and there is no friction, this graph will trace a circle as the mass of fluid settles 
into a stable path of oscillation. 

The purpose of creating the two-dimensional model is to investigate the mathemat¬ 
ics and physics of the shallow water equations before moving on to modeling the three- 
dimensional shallow water equations. Furthermore, because of the simplicity of the two- 
dimensional shallow water equations, this first model incorporates more features than will 
be implemented in the three-dimensional model within the scope of this project. Factors 
such as the residual velocities at each depth layer, wind friction, and variable viscosities have 
not been included in the three-dimensional model. 


8.1 Operation of the Two-Dimensional Model 

The two-dimensional model is implemented as a MATLAB function that can be run 
from the MATLAB base workspace without any input arguments. The initial conditions are 
specified through the use of another function call that is designed to be user-modifiable. The 
actual calculations are carried out within a loop in the main function of the model. 


The most important single line of this function is derived directly from (75), which 
is the first step of the two-step operator splitting method established in (76). This first-step 
differential equation (75) is 


M 



-F pres - F pot + Mu + M w (u w - u) - M fr (u). 


Since the two-dimensional model solves this equation using an implicit scheme, the derivative 
of u is approximated by the backwards Euler method. Using this approximation, (75) 
becomes 

{ — W n \ — 

M f-—-j = —F pres — F pot + (M — M w — Mf r ) u" +1 + M w u w 

where superscripts represent the time step at which a variable is evaluated. Solving this 
equation for u n+1 yields 


M — dt (M - M w 


M fr ) 


u 


71+1 


= dt 


pres 


F P ot + M w u w + Mu r 
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This equation is used to compute the velocity of each tracked particle in a certain time 
step using the velocity of the previous time step. The various matrices for mass, viscosity, 
and friction coefficients must be created within each time step previous to evaluating this 
equation. 

When implemented in MATLAB, the value of u n+1 is determined using linear algebra; 
specifically, MATLAB’s backsolve feature. The code that solves this equation is 

U = ((MM - MMVV + dt*(MMWW + MMFF)) \ (-dt*POT' - dt*PRES' + ... 

dt * MMWW * UW ' + MM*U') ) '; . 

In this code, U represents the vector of velocities for all tracked particles. It is updated with 
new values when this line of code finishes execution. The matrices M : M, M w , and Mf r 
correspond to MM, MMW, MMWW, and MMFF, respectively. The force vectors F pres and F pot 
become PRES and POT. The one discrepancy between the equation and the implementation 
is that the viscosity matrix-M and MMVV-is multiplied by dt in the equation but not in the 
implementation. This is because M as defined in (73) is already multiplied by a factor of 
dt. Rather than multiplying by this factor and then later dividing it out, it is more efficient 
to simply omit this factor in the code. 

After solving for new values of U in each time step, we determine the vector of each 
tracked particle’s position, represented in the code by X. This is the second step of (76). 
Recall that the differential equation we solve in this step is 

dx 

— = u. 
dt 

This is implemented in the model as 
X = X + U * dt; 

After this step is complete, the last piece of information of interest to the model is 
the water column. When the method for changing from a system of continuum equations to 
discrete equations was presented in Section 7, the relationship between water column and 
mass within each fluid element was established in (56) to be 


— 

%i- 1-1 

This relationship is established at time zero at which time the mass, m°, within each fluid 
element is fixed. For all other times this relationship still holds; h % and Xi + 1 — Xi change with 
each time step while m, stays constant. Therefore, in the implementation of the model, the 
water column is calculated by the code 


H = M./diff(X) ; . 
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The diff function in MATLAB computes the difference between adjacent elements of the 
vector X. For a vector holding the elements X\ to x n+ 1 , diff (X) returns a length n vector 
where each element has the value Xi + 1 — x t . 

The three lines of code in this chapter are the heart of the implementation of the two- 
dimensional model. They lie within a “time-stepping loop” that counts up time steps from 
time zero when the program is started. The rest of the code within the loop is devoted to 
setting up the various matrices and forces needed to compute the particle velocities and then 
plotting the resulting data. Matrices M and M are constant and are set at the beginning 
of the program before it enters the time-stepping loop. The pressure and potential terms 
are set in sections of the code labeled after them, while M w and Mf r are set in the wind 
and friction section. Within the loop two indices are maintained: t, the time index which 
represents model time, and i, the counting index which is the number of iterations through 
which the loop has passed. Each time the counting index advances by 1, the time index 
advances by dt. 

Visualizing the data has no real part to play in the actual operation of the two- 
dimensional model. The model will create numerical output representing the velocity and 
position of fluid particles and water column of fluid elements from any valid input and 
without the help of any plotting. However, visualization provides the human user with a 
way to immediately understand the data that is being generated from the model without 
having to read and compare huge vectors of raw numbers. As mentioned in the previous 
section, there are two graphical outputs from the two-dimensional model: the water column 
versus x plot and the center of mass versus average momentum plot. Plotting does not have 
to be carried out by the program during every single time step. Instead plotting is done at 
a specific time interval, which can be specified by the initialization file. 

The initialization takes place before the program enters the time-stepping loop. A 
separate function is run that creates initial conditions based on user specifications. Through 
its development process, the initialization program has expanded to encompass more and 
more user-input so that now almost every aspect of the program is modifiable, from the 
magnitude of viscosity and friction, to the initial shape of the bathymetry and wetted area, 
to the manner in which functions are defined as inputs. The initialization file is run at the 
beginning of the model’s executable and it loads its state date into a structure which is then 
passed to the model itself; the model updates this data structure as it steps through time. 
Once the model is done running, it compiles all the data back into the state structure which 
it returns as function output. This data structure can thus be used as input for another run 
of the model. 

The complete code for the two-dimensional model is included in (C) and the initial¬ 
ization code is included in (B). 
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8.2 Conclusions from the Two-Dimensional Model 

The two-dimensional model is a stepping stone that we created to gain better un¬ 
derstanding of the mathematics needed for the three-dimensional model. Because the two- 
dimensional model has fewer governing equations it is simpler to trace through the series of 
transformations needed to bring this model to implementation. We experimented with our 
novel ideas on this simpler model before tackling the three-dimensional problem that is the 
ultimate goal of this Trident project. 

For first few steps of transformation of the two-dimensional model’s governing equa¬ 
tions, we followed the precedent of the Princeton Ocean Model. The non-dimensionalization 
and s-Coordinate transformation both originated with POM. However, we successfully in¬ 
troduced several novel ideas in our derivation. One new concept is the way we handle the 
residuals after depth-averaging the conservation of mass and momentum equations in section 
6. We replace the evolution equation by the local equilibrium solution, thereby removing a 
degree of freedom from our model. This will prove to be particularly useful for the more 
complicated three-dimensional model. 

Another innovation is using a Lagrangian particle tracking approach for our imple¬ 
mentation rather than the Eulerian method used in POM. By solving the problem on a 
domain with a boundary fixed by the locations of the fluid-marker points that we track, we 
obtain information about the free boundary. This makes our model ideal for determining 
if a location is wetted or not and also helps to make our model more accurate for shallow 
water estuaries and enclosed bays. 
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9 Moving from the Two Dimensions to Three Dimen¬ 
sions 

We now move from the shallow water equations with one lateral dimension to those 
with two lateral dimensions. A model of these equations has already been created by J. M. 
Greenberg, an adviser on this project. The objective of this part of the project is to take 
this model—which is implemented in serial for one processor—and alter it so that it can 
run in parallel on multiple processors. To do this, the serial code must first be modified to 
make it amenable to efficient parallelization. This entails both dividing up the processor 
time required for the calculations in the model so that each processor has a roughly equal 
load and providing a means for communicating data between processors. 

The mathematics of this model are similar to those of the two-dimensional equations 
from sections 3 through 7, except that there are now two lateral coordinates, x and y , in 
addition to the vertical coordinate, z. The three dimensional versions of the basic physical 
equations of (6)-(8) are 


u x + v y + w z = 0, (77) 

u t + uu x + vu y + wu z +p x = {<T\\) X + ( 0"12 ) y + (offis), + f v , (78) 

v t + uv x + vv y + wv z +p y = (a ii)^ + {(J 22 ) y + ( 043)2 - fu, (79) 

w t + uw x + vw y + ww z +p z = (< 731)3 + (o- 23 ), y + (a 33 ) z - g, (80) 


where density is a constant equal to one and v is the velocity in the direction of the new 
lateral coordinate y. In addition, the Coriolis acceleration is now an important factor; fv and 
fu represent its effect in this model. From this equation we proceed to a set of differential 
equations. Instead of having just one equation where we solve for % as a function of x and 
u , we obtain two equations; one for and one for both in terms of x, y, u, and v. 

The one important difference between the implementation of the two- and three- 
dimensional models is the “geometry”—the decomposition of the lateral space in the fluid 
domain into particle locations and fluid elements. This difference is discussed in greater 
detail in the next section. Another significant difference is that the three-dimensional model 
abandons the implicit Euler method for solving the differential equations in favor of the faster 
explicit Euler method. Because the number of particles tracked by the three-dimensional 
model is generally much larger than in the two-dimensional model, speed is of greater im¬ 
portance. In addition, this project’s implementation of the three-dimensional model ignores 
some complicating factors incorporated into the two-dimensional model. This removes one 
of the chief advantages of the implicit method in the two-dimensional problem; namely, that 
the implicit method is more stable with complex or unusual initial conditions. 


While the equations are being solved by a different method, the process for achieving 
a solution is still similar. The three-dimensional model still uses the same two-step operator 
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splitting method as the two-dimensional model. In the first part of each time step x and y 
are held constant and u and v are updated; in the second part u and v are held constant 
and x and y are updated. The three-dimensional model uses the same time grid as the 
two-dimensional; in both versions dt is the length of the time interval between time t n and 
t n+ 1. The two steps in the three-dimensional model can be described by 


1VT — = TU 


Step 1 

Step 2 

t n <t< t n+1 

t n <t< t n+1 

dx. dy q 

dt dt 

^ = u ^ = V 
dt ’ dt v 

u,v), M y ^ = <7(x,y,u,v) 

M Ai — M — = 
iVlx dt. iV1 y dt 


The mass matrices M x and M y differ from the mass matrix M of the two-dimensional model. 
The two-dimensional model used tri-diagonal matrices and required linear algebra to solve. 
Because we are trying to avoid the use of linear algebra, M x and M y are appropriately 
chosen diagonal matrices in the three-dimensional model. 



Figure 4: One-dimensional grid 


9.1 The Geometry of the Three-Dimensional Model 

A Lagrangian solution to the shallow water equations requires a grid where particles 
are identified at time zero and then tracked through all subsequent time steps. While the 
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two-dimensional model has a simple linear geometry, the three-dimensional model is much 
more complex. Since there is only one lateral degree of freedom in the two-dimensional 
model, the “grid” was simply a line with particle locations that could move back and forth 
across the grid in the one degree of freedom, as seen in (Figure 4). The n + 1 particles divide 
this one-dimensional grid into n fluid elements. 



Figure 5: Two-dimensional grid 


For the three-dimensional model, we use the two-dimensional grid depicted in (Figure 
5). In this grid, there are two types of nodes which are both analogous to the particle 
locations of the two-dimensional model. Some nodes are connected to their four neighboring 
nodes in the four cardinal directions; these nodes are designated by circles in the figure. 
The other type of node is connected to eight neighbors-the four cardinal neighbors and four 
more diagonal neighbors-and is represented by a square. We cannot simply make every node 
an eight-node because then there would be additional points where edges intersected that 
were not in the initial node grid; specifically, where the diagonal edges of adjacent nodes 
crossed each other. The four- and eight-node grid allows for a greater number of edges while 
demanding that edges meet only at the nodes. 

The significance of the edges is in the geometry of the triangles that they form in the 
grid. Each interior node of the grid has a square split into two triangles to its upper left. 
These two triangles are identified by the node to their lower right; the nodes are numbered 
sequentially as in (Figure 6). The triangles are divided into two classes depending on their 
spatial orientation, numbered either “1” or “2”. Triangle numbers stay the same over vertical 
and horizontal edges lines, while alternating over diagonal edges. This is shown in (Figure 
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7). We can see from this figure that each square has both a 1-type triangle and a 2-type 
triangle in it. This allows us to distinguish the two triangles from each other. Each triangle 
is identified by the node that is below it and on its left, as well as by its type number. 



Figure 8: Two-dimensional grid with one pair of triangles highlighted 


With this numbering scheme for triangles, we now have a way to label the edges of 
each triangle. For example, in (Figure 8), the darkened triangle can be identified as the 
“1-type” triangle belonging to node (3,4). When this model is implemented, it will be 
important to keep track of the changes in x, y , u, and v along each of the edges in the 
grid; this trianglular numbering scheme allows us to have a unique identifier for each edge. 
Another important piece of information to store is the area of each triangle. Because this 
is a Lagrangian scheme, the mass in each triangle at time zero will be constant for all time. 
In the two-dimensional model, this was also the case, and the water column in each fluid 
element was calculated by dividing the constant mass in that element—as established in the 
initial conditions at time zero—by the distance between the two fluid particles that form the 
boundary of that fluid element. In the three-dimensional version, the water column in each 
triangular element will be a constant that can be calculated from the constant mass in the 
triangle and the area of the triangle as it deforms through time. 

The grid depicted in (Figure 5) has an eight-node—that is, a node that is connected 
to eight neighbors—at all four corners. If the nodes are numbered in a m+1 x n+1 square, 
then both m + 1 and n+1 must be odd in order for all corners to have an eight-node. The 
grid is implemented this way in order to be consistent, such that there will always be a set 
triangle geometry on the edges. By demanding a set geometry, the implementation of the 
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model becomes much easier: there is no longer any need to test the types of node at any 
given point. Instead we can be sure that all nodes where both m +1 and n +1 are both either 
odd or even are eight-nodes while all nodes where one of m+1 and n +1 is even and the other 
is odd are four-nodes. For example, (1,1) and (4,2) are both eight-nodes while (1,2) and 
(4, 3) are four-nodes. It is important to know whether a node is an eight-node or four-node 
because that tells us about the orientation of the two triangles that belong to that node. 
For an eight-node, the hypotenuse of the two triangles has a positive slope, that is it goes 
from the bottom left of the square to the upper right. The hypotenuse of the four-node has 
a negative slope. 

9.2 Implementing the Three-Dimensional Model in Serial 

The serial version of the two-dimensional model solves the final form of the shallow 
water equations explicitly. This is done because the linear algebra needed to use the implicit 
Euler method becomes extremely time-consuming as the number of edges connecting nodes 
increases with the two-dimensional grid. Despite using the faster explicit Euler method 
to solve the differential equations of the operator splitting method, the three-dimensional 
model still runs significantly slower than the two-dimensional version. One reason is that 
more processor time is spent calculating the differences between nodes along the various 
edges; this data is needed to calculate pressure and viscosity terms. In the two-dimensional 
model there is a maximum of two edges per node along which differences must be calculated 
for the variables of interest: depth-averaged velocity, position, and water column. In the 
three-dimensional model, as noted in the previous section, there are up to eight edges per 
node; this is what causes the linear algebra to be impractical with the implicit Euler method. 
In addition, there are more variables of interest: depth velocity is represented by both a u 
and v component and position is represented by both x and y. The calculation of differences 
in water column is even more complex. In the two-dimensional version, the fluid elements 
are connected only at tracked particle locations, so it is only at these locations that we are 
interested in finding the difference in water columns. In the three-dimensional version the 
fluid elements are the triangles pictured in (Figure 5). There are approximately two triangles 
per node, and each triangle has three edges that have a unique water column difference. 

Because the three-dimensional model uses the explicit method to solve its differential 
equations, care must be taken to ensure the stability of the model. In the two-dimensional 
model, the implicit method for numerical approximation was chosen to help ensure stability. 
Since the three-dimensional model uses the explicit method, it is necessary to choose time 
steps that are sufficiently small to ensure that there is no “cross-over.” Cross-over is an 
error that occurs when the nodes of (Figure 5) move in such a way that some of the triangles 
overlap each other. The equivalent error in the two-dimensional model occurs when a particle 
that started with a x value greater than an adjacent particle gets a position value that is less 
than that adjacent node. Carefully choosing time-steps that are sufficiently small can prevent 
this error from occuring: in practice stability has been ensured when the ratio between dtt 
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and the dx and dy of the initial grid is 0.02 or less. 

The increased complexity of the three-dimensional model in serial causes greatly 
increased runtime; this is the motivation for the parallelization efforts of this project. The 
computationally intensive parts of the model must be organized in a way that makes it 
easier to divide the workload of the three-dimensional model over multiple processors. The 
serial implementation of the three-dimensional model can be broadly divided into three 
parts. Recall from the two-dimensional implementation that there are three lines of code 
which constitute the “core” of the model; these lines compute the depth-averaged velocity, 
position, and water column and exist within the time-stepping loop. The other code within 
the loop computes the various forces and matrices that are necessary to evaluate the three 
core lines. The three-dimensional code is constructed similarly. There are a few core lines 
that compute the three lines of interest, and these lines are inside the time-stepping loop 
along with code that calculate the necessary forces and matrices-which are much more 
complex in the three-dimensional model as mentioned above. 

The three parts into which the three-dimensional model can be divided are everything 
inside of the time-stepping loop, the loop itself and everything outside the loop. The calcu¬ 
lations inside the time stepping loop can be treated as a single-step calculation that takes 
the three variables of interest on a two-dimensional grid and returns these variables on that 
same grid one time step later. This can be treated as a separate function, and implementing 
the serial version in this way greatly eases the transition to the parallel version. Note that 
the grid input into this one-step function must follow all the rules established for the grid 
in the last section; specifically, that it must be rectangular and have an odd length in both 
dimensions in order to ensure that there is an eight-edge node on all four corners. 

The loop itself is almost trivial in the serial version of the three-dimensional model. It 
only needs to incorporate the same two counting indices from the two-dimensional model —t 
and i —along with the logic for advancing the loop and ending it at the end of the runtime. 
The last part is the setup code that runs before entering the loop. This can also be treated 
as a separate self contained function. 


9.3 Implementing the Three-Dimensional Model in Parallel 

To produce efficient parallel code, work must be divided equally among several processors. 
A simple way to do this is to divide the two-dimensional particle grid evenly between the 
number of available processors. This then presents the problem of communication: each 
node has a separate set of differential equations to solve-specihcally the solutions to 

, „ du _ 

M x — = .F(x,y,u,v) 


My 


dv 

dt 


and 


£( x ,y,u,v) 
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which constitute the first step of the two-step operator splitting paradigm that this model 
uses. Both of these equations use vector variables to represent a set of n + 1 differential 
equations, each equation operating at one node on the two-dimensional grid. A sub-set 
of these equations and their associated nodes will be assigned to each separate processor. 
However, to solve any of the differential equations at a node, it is necessary to have data from 
all nodes connected to that node. If the node in question is on the edge of an area that has 
been split from the whole and assigned to a certain processor, then it may need data from 
nodes that have been assigned to a different subset on a different processor. Overcoming this 
problem requires a scheme for passing data between parallel processors, and this requirement 
drives the way in which the sub-domains are divided from the whole two-dimensional grid 
and assigned to processors. 

Communication is, in general, more time consuming than processing. To optimize the 
parallel code, communication must be avoided as much as possible; therefore, the process for 
splitting the two-dimensional grid is designed to keep the amount of communication necessary 
small. Furthermore, each parallel processor should be assigned an equal share of the work 
to perform. Each time a processor pauses to communicate, it must wait for all processors to 
complete the computational tasks they had been assigned. Assuming that each processor is 
of roughly equal speed—an assumption that holds true for the Naval Academy’s cluster— 
each processor should be assigned an equal sized sub-domain of the whole two-dimensional 
grid and equally sized sub-set of the differential equations. 



Figure 9: Two-dimensional grid divided into square sub-domains 


Two of the easiest ways to divide the two-dimensional grid into sub-domains are by 
squares as depicted in (Figure 9) or by rectangles as in (Figure 10). If the goal is simply to 
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Figure 10: Two-Dimensional grid divided into rectangular sub-domains 


divide in a way that results in minimal communication, the square is the superior method. 
For the rectangular method on an m+1 x n+l grid with P processors, there are n +1 (P — 1) 
edges that cross sub-domain edges. For the same m+1 x n+l grid with sub-domains assigned 
to P processors by the square method there are (m + n + 2) (^/P — 1^. Despite requiring 
more communication, the rectangular method has several advantages. Firstly, it is divided 
in only one coordinate. Thus, only the n columns in the x —or the y direction depending 
on the orientation of the grid—need be divided among the processors. This allows each 
processor to evaluate over the 1 : m in the non-divided direction. This not only simplifies 
the code for solving differential equations at each node, but also greatly simplifies the code 
for communication between processors. 

A second advantage of the rectangular method is in the pattern of communication. 
In a cluster using the square method for creating sub-domains, each processor may need to 
communicate with as many as eight other processors. With the rectangular method, each 
processor is guaranteed to only have two communication partners at most. Depending on 
the configuration of the hardware of the cluster, this can be a great advantage. For clusters 
organized in a lattice or a ring, the necessity to only communicate with two adjacent proces¬ 
sors eliminates the need for any communication to pass through multiple processors in order 
to get from its source to its destination. Despite the fact that the Naval Academy’s Beowulf 
cluster uses an InfiniBand Ethernet connection that allows all nodes to communicate with 
each other equally, this project’s three-dimensional model uses the rectangular method to 
create sub-domains because of its increased viability on other types of clusters and espe¬ 
cially because it greatly simplifies the process of writing code for communication between 
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processors. The rectangular shape of each sub-domain is the first constraint in the scheme 
for splitting the whole two-dimensional grid. 

Determining that this project will use the rectangular method to create sub-domains 
from the whole two-dimensional grid sets the shape of each sub-domain. The next step 
is to determine the size of these sub-domains. The driving factor when determining the 
size is the three parts into which the serial implementation of the three-dimensional model 
was divided. Recall that all of the calculations—including the core code that computes the 
three variables of interest—within the time-stepping loop can be organized into a discrete 
single-step function. The parallel version of the code will use this same function completely 
unchanged from the serial version. In order to support this, the parallel version of the code 
must be able to pass this function data from the nodes of a two-dimensional grid that has all 
the same properties as the whole grid. The rectangular shape of each sub-domain already 
ensures that the m+1 dimension of each sub-domain is the same as in the whole m+1 x n+1 
grid. The n+l dimension of each sub-domain must have an odd number of nodes so that the 
four corner nodes of each sub-domain are eight-edge nodes. This is the second constraint for 
the sub-domains. 
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Figure 11: Sub-domain scheme with adjacent columns as used in three-dimensional model 


A third constraint on each sub-domain is that each sub-domain should contain an 
extra column of nodes that “belong” to each of the up to two adjacent sub-domains. As 
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in (Figure 11), the left-most sub-domain should contain at least one column of nodes from 
the second left-most sub-domain. The reason for this is that in order for each node that 
“belongs” to a processor to be correctly updated, that processor must have access to the 
data from all adjacent nodes. These two constraints are enough to develop a scheme for 
splitting sub-domains from the whole two-dimensional grid. 

Let the number of parallel processors that are available to the model be P = 2 P 
and let the number of nodes along the m + 1 dimension of the two-dimensional grid be 
m + 1 = 2 m + 1. While these two assumptions hold, each processor will be responsible for a 
sub-domain of 

nM 

_ _ c\M—p 

2 p ~ 

columns of n +1 nodes, except for the right most sub-domain which will “own” the one extra 
column for a total of 2 M ~ P + 1 columns of n + 1 nodes. This division is pictured in (Figure 

10) . To fulfill the second constraint each sub-domain will have attached to it a column of 
nodes from each adjacent sub-domain. To fulfill the third constraint, each sub-domain will 
have one extra column from the adjacent sub-domain to the left appended to it as in (Figure 

11) . As can be seen in the figure, each of the sub-domains now has a rectangular shape, has 
eight-edge nodes at all four corners, and contains all the data necessary for the processor to 
perform single-step calculations on all the nodes that processor owns. Each processor has 
been assigned 2 M ~ P + 3 nodes except for the Erst processor which only has 2 M ~ P + 1. It 
would be better for each processor to be assigned the exact same number of nodes, thereby 
allowing each processor to run exactly the same single-step computation function. Therefore, 
the first processor is assigned two extra columns of nodes to the right of it, bringing its node 
count even with all the other processors. 

The parallel version of the model first runs a modification of the set-up code that 
began the serial version. It then moves into the loop. During each iteration of the loop, it 
runs the same single-step function that was run in the serial version; making this possible is 
the prime motivation for the scheme for grid splitting. The last piece that must be added to 
make the parallel version work is a section within the loop that communicates the necessary 
data between processors. The geometry establishes that each processor needs to be passed 
only two columns of node data—or one column for the two end processors—to run each time 
step. Therefore all interior processors must pass two columns of node data and receive two 
columns of this data, while the end processors must pass and receive one column. 

The extra columns of nodes that were appended to each sub-domain to comply with 
geometry requirements do not need to be passed. There were initially concerns that because 
this data is not updated to reflect the changes taking place in other parts of the model, that 
it might eventually become so divergent from the data that is passed that it would cause 
errors. The most likely error that would be caused would be the “zero area triangle” error in 
which one of the triangle’s area approaches zero and therefore its water column approaches 
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infinity. However, it has been determined experimentally that although this data will not 
be updated accurately, it will remain “reasonable” enough to not cause any errors in the 
program; no run of the three-dimensional model has ever caused this error. 
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10 Results of Test Data and Conclusions 


Table 1: Results of Testing 


Domain Size 

Serial 

4 Proc. 

Rect. 

4 Processor 

Parallel 

8 Proc. 

Rect. 

8 Processor 

Parallel 

16 Proc. 

Rect. 

16 Processor 

Parallel 

65 x 65 

0.7953 

0.2839 

0.3766 

0.1939 

0.2976 

0.1461 

0.3158 

129 x 129 

4.1803 

0.8650 

0.9338 

0.4878 

0.5743 

0.3205 

0.5210 

257 x 257 

24.5459 

4.7274 

5.0877 

1.8032 

1.9852 

0.9316 

1.1136 


The data in (Table 1) was generated through test runs on the Naval Academy’s 
Beowulf cluster with the MATLAB Distributed Computing Toolbox. The tests are run on 
square two-dimensional grids whose sides are of length 2 n + 1. Three grid sizes were tested: 
n — 6, n — 7, and n = 8. The code for both the serial and parallel versions were modified 
from the versions attached in (D) and (E) to track runtime. Note that the times in (Table 
1) are not processor times but “real” times. The code starts tracking time for each test 
after the initialization of data is done and right before the time-stepping loop begins. In the 
parallel version, this means that the initial grid and data are distributed to all processors 
before the timing starts. 

The domain size is the number of nodes along each edge of a square grid created for 
the test. The column of data labeled “Serial” is that runtime for that domain size in the 
serial version. The data labeled “Parallel” is the runtime for P parallel processors. The data 
labeled “Rect.” represents the runtime for one single processor when the domain is divided 
between P parallel processors. This represents the theoretical maximum speed of the model 
without any communication. The data in (Table 2) is calculated from the data in (Table 1) 
and shows the speedup factor for each domain size and number of processors, as well as the 
percentage of time spent communicating. 


Table 2: Speedup and Communication Time 


Domain Size 

4 Processor 

Speedup 

4 Processor 

Comm. Time 

8 Processor 

Speedup 

8 Processor 

Comm. Time 

16 Processor 

Speedup 

16 Processor 

Comm. Time 

65 x 65 

2.11 

24.6% 

2.67 

34.8% 

2.52 

53.7% 

129 x 129 

4.48 

7.4% 

7.28 

15.1% 

8.02 

38.5% 

257 x 257 

4.82 

7.1% 

12.36 

9.2% 

22.04 

16.4% 


The results of the test data show that the implementation of the three-dimensional 
model in parallel is successful. For a sufficiently large grid, increasing the number of pro¬ 
cessors available will increase the speed of the model. There is a reduction in efficiency as 
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the number of processors increases. This is result of increased communication between the 
processors. However, this is partially counterbalanced by a significant reduction in runtime 
for smaller domains. The ratio of runtime between larger and smaller domains is greater 
than the ratio between the number of nodes in each domain. This means that while one 
domain may be half the size of another, the processing time for the smaller domain will 
be less than half of the time for the larger domain. We suspect that this is because larger 
domains require more memory usage, and make proportionately less use of caching when 
performing calculations. 

We see that the speedup on larger domains can exceed the number of processors. This 
can be seen in the tests on the grid of size n = 8, where the speedup with 8 processors is over 
12, and the speedup with 16 processors is over 22. However, once the number of processors 
becomes too large, the losses introduced by communication overcome the advantages gained 
by using smaller domains. In the grid of size n = 6, 16 processors actually run slower than 8 
processors because communication time takes up such a large percentage of the total runtime. 

The results of this test data confirm that the methodology for splitting work between 
processors used in this project is a success. The deliverable three-dimensional model created 
by this Trident project is a proof-of-concept for the domain-splitting and communication 
strategy we introduce. As the complexity is added to the three-dimensional model we have 
created, the strategy for parallelization should remain valid. 
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A Surface Traction Derivation 


The normal to the surface z = a(x) is 

1 


n s = 


—A 


y/l + A 2 a| L 1 
and the normal to the surface A = a(x) + h(x,t) is 

1 

^a+h 


(A-l) 


A (a,x + h s ) 
1 


-y/l + A 2 (a £ + hs) 

Using the values of a n and a 13 from (17) and (A-l) and (A-2) we find that 

- 2 ea s Ux(x, a,t) + e (wx(x, a,t) + ^u 2 (x, a, t)) 


(A-2) 


(011) ^13) ‘ n a ~ 


and 


(0ii,013) • n s+/l — 


y/l + A 2 a| 

2e ( a,x + hx ) a + d, f) + e (%(!, a + h, t) + ^^(x, a + d, f)) 
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^/l + A 2 (a,x + /i^) 


(A-4) 


Let J be the physical Cartesian stress tensor 


or in non-dimensionalized form 
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(A-5) 


The traction vectors are 

and 


T a = J n s 


^a+h J ' ^a+/i' 


We evaluate the traction vectors using (A-l), (A-2), and (A-5) to obtain 


and 
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Note that the first components of T a and T (1+/ - are (dn, di 3 ) • n a from (A-3) and (dn, d 13 ) • 
n a+a from (A-4), respectively. The tangential components of these traction vectors are found 
by the dot product of the traction vectors and a unit vector tangential to the surface as in 

$ 5 , = (j ' n a) ■ ‘S’a+h = (J ' n a+/i) ' ^a+tr 

The tangential components of the traction vectors (A-6) and (A-7) are 

S ' a = l (1 + A 2 d 2 ) [j ^ " A ^ + A “ AV )] (A " 8) 

and 

S a +h — —/ - — 7 — -=—v [% (l — A 2 (b 2 + h 2 )) + \w s (l — A 2 (b 2 + h 2 ))1 (A-9) 

a+h l (1 + A 2 (d 2 + h 2 )) La v v n xK v n \ K ’ 
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B MATLAB Code: initialize2d.m 


function struc = initialize2d() 

%***************************USER INPUT AREA**************************** 
%NOTE: ALL USER INPUT FUNCTIONS MUST BE DESIGNED FOR MATRIX OPERATIONS. 
%THUS INSTEAD OF *, /, AND ", THE INPUT SHOULD BE IN THE FORM ./, 

%AND . " UNLESS THE OPERATION IS BETWEEN A VARIABLE AND A SCALAR. 

%Settings for this simulation. A number of features can be turned on or 
%off. THese features are set below, with either 'true' for on, or 
%'false' for off. 

%Timing - If set to 'true', this run will measure the execution time 
%for the number of time steps specified in 'timesteps'. If 'false,' 
%timesteps is ignored and runtime is set by runtime variable below, 
timing = false; 
timesteps = 100; 

%Movie - If set to true, this run will be recorded as a movie. 

%name is the name of the file to which the movie will be stored 

mov = false; 

name = 'moviel.avi'; 

%Residual - If set to 'false' only the depth averaged mean will be 
%used in calculations, and none of the variation at different depths 
%will be considered 
resid = true; 

%Friction - If set to 'false' then there will be no friction terms in 
%the calculations. Note that setting friction to 'false' automatically 
%sets the residual to 'false' as well. Also note that the upper and 
%lower friction constants can be turned off individually with the Kw 
%and Kf constants below, 
friction = true; 

% Set N or dx. if isN is true, then val is N, if isN is false, then val 
% is dx. 
isN = true; 
val = 100; 

% Set the function defining the bathymetry as a function of x 

% If bathymetry is an inline function, bath_type is true, if it is an 

% m-file function, then bath_type is false. 

bath_type = true; 

bathym = '(x."2)/2'; 

% Set the domain of the bathymetry. This domain should be the absolute 
% boundary of the problem. For now the only purpose of this is to set 
% axes for plotting. 
hi_b = 3.5; 
lo_b = -3.5; 

% Set the initial domain of the fluid 
hi_f = 2; 
lo_f = 1; 

% Set the function defining the initial water column as a function of x 
% If water column is an inline function, water_type is true, if it is 
% an m-file function, then water_type is false. 
water_type = true; 
water = '0*x+5/3'; 
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% Set the function defining the initial velocity of the fluid as a 
% function of x 

% If initial velocity is an inline function, vel_type is true, if it is 
% an m-file function, then vel_type is false. 
vel_type = true; 
vel = ' 0*x+0'; 

% Set constants for surface frictions. Kw is the upper surface friction 
% and Kw is the lower surface friction. If these terms are set to zero, 
% then there will be no effect from that friction source. 

Kf = 0.002; 

Kw = 0.0005; 

% Set wind as a function of x and t FOR NOW JUST X 

% If wind is an inline function, wind_type is true, if it is a m-file 
% function, then wind_type is false. 
wind_type = true; 
wind = ' 0 *x-l'; 

% Set the value of lambda 
lambda = .001; 

% Set the value of mu. Mu is the percentage between 0 and 1 of the 
% damping to be used 
mu = .8; 

% Set the ratio dt/dx. A typical value is 0.1 
rat = 0.1; 

% Set the runtime for this model, 
time = 100; 


%************************ERROR CHECKING AREA*************************** 


%*************************CALCULATION AREA***************************** 

% Set the features to be turned on and off 
struc.res = resid; 
struc.fric = friction; 

% Set the bathymetry and its domain, 
if bath_type 

struc.bath = inline(bathym,'x'); 

else 

struc.bath = bathym; 

end 

struc.hi = hi_b; 
struc.lo = lo_b; 

% For a given N, calculate the x-array: constant mass, variable dx 
if isN == true 

struc.n = val+1; 
if water_type 

wat = inline(water,'x') ; 

else 

wat = water; 

end 

el_mass = quad(wat, lo_f, hi_f, l.e-10) / val; 
int_dom = lo_f: (hi_f-lo_f)/10000:hi_f; 

= cumtrapz(int_dom, feval(wat, int_dom)); 


cum 
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struc.X(l) = lo_f; 
j = 1; 

for i = 1:10000 

if cum(i) > j * el_mass; 

j = j + 1; 

struc.X(end+1) = int_dom(i); 


end 


end 

struc.X(end+1) = hi_f; 


approx_dx = (hi_f - lo_f) / val; 
struc.dt = rat * approx_dx; 


% For a given dx, calculate the x-array: constant dx 
else 

% Since not every dx will divide evenly into the domain, round off 

% N to the next highest whole number and recalculate dx. 

struc.n = ceiling((hi_f-lo_f) / val) 

struc.dx = (hi_f-lo_f) / (struc.N-l) 

struc.dt = rat * struc.dx; 

struc.X = lo_f:struc.dx:hi_f; 

end 

% Set the initial h vector. First evaluate 'water' at all points in the 
% x-array, then find the midpoints between these values, leaving N-l 
% datapoints in h. 
for i = 1:length(struc.X)-1 

struc.M(i) = quad(wat, struc.X(i), struc.X(i+1), l.e-10); 

end 

% Set the initial u vector. Evaluate 'vel' at all points in the 
% x-array. 
if vel_type 

struc.U = feval(inline(vel,'x'), struc.X); 

else 

struc.U = feval(vel, struc.X); 

end 

% Set friction constants of the model. 
struc.k_lo = Kf; 
struc.k_hi = Kw; 

% Set the wind function for the model 
if wind_type 

struc.uw = inline(wind,'x'); 

else 

struc.uw = wind; 

end 

% Set the lambda value for the model 
struc.lam = lambda; 

% Set the mu value for the model 
struc.moo = mu; 

% Set the runtime of the model, 
if timing 

struc.time = timesteps; 
struc.timing = true; 

else 

strue.time = time; 
struc.timing = false; 


end 



%Determing if movie or not 
struc.mov = mov; 
if mov 

struc.name = name; 

else 

struc.name = 

end 
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C MATLAB Code: shallow2d.m 

function out = shallow2d() 

%close all other windows 
close 

%%Autoruns initialize.m to get setup data. Alternate option is to allow user 
%to specify input file in function definition, 
data = initialize2d(); 

%used to time execution of the program 
if data.timing 
tic; 

end 

%Get initialization data from initialize.m 

resid = data.res; 

friction = data.fric; 

bath = data.bath; 

n = data.n; 

dt = data.dt; 

X = data.X; 

M = data.M; 

U = data.U; 
hi_b = data.hi; 
lo_b = data.lo; 
kw = data.k_hi; 
kf = data.k_lo; 
wind = data.uw; 
lambda = data.lam; 
mu = data.moo; 
runtime = data.time; 

%Set a vector of the mass in each interval and total mass in problem 
H = M. /diff(X); 

%%initialize time and counting indices to zero 
t = 0; 
i = 0; 

% sundry variables set to zero: HH and XX for formatted display, xc and vc 
% are centers of mass and momentum 
XX = 0; 

HH = 0; 

XC = 0; 

VC = 0; 

T= .1/dt; 
mov = 0; 

% constants for plotting 
%axes for mass vs. momentum 
AX = - 2 : .1:2; 

AY = AX*0; 

%x and y for plotting potential function 
X_POT = data.lo:.1:data.hi; 

Y_POT = feval(bath, X_POT); 


%% calculate (time independant) terms 
% Mass matrix 

MM = spdiags([[M 0]'/8,3* ( [0 M] + [M 0])'/8,[0 M]'/8],-l:l,n,n); 
% Viscosity matrix 
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MMV = spdiags ( [ [M 0]'/8,([0 M] + [M 0] ) ' /—8, [0 M]'/8],-l:l,n,n); 

% Total mass 

mass = sum(MM * ones(size(1:n))'); 

%% Looping 

%If timing then count by index numbers, otherwise count by elapsed time 
if data.timing 

runtime = runtime * dt; 

end 

while t < runtime 

% update time increment to t_n+l and increment counting index 
t = t + dt; 
i = i + 1; 


% WE UPDATE U TERM BY TERM 

%%MASS AND VISCOCITY MATRICES 

%MM has already been set 

MMVV = MMV * mu; 

%END MASS MATRIX 

%%PRESSURE TERMS 

if resid == true 

%This section evaluates the integral of the square of the residual at 
%each point x. This becomes one of the pressure terms in the solution 
%for the depth averaged mean velocity. 

%These next calculations must be evaluated at x_i-l/2 and x_i+l/2 for each 
%point x_i. This means there will be n+1 points of evaluation where 
%n is the number of x points. For the two endpoints at x_l/2 and 
%x_n+l/2, we will estimate their value by x_l and x_n, 

%respsectively. 

U_AVG = ([U(1) U] + [U U(end)])/2; 

%We set X_AVG to equal all the midpoints so that we can find the 
%wind function for those points. 

X_AVG = ([X(l) X] + [X X(end)])/2; 

UW_AVG = feval(wind, X_AVG); 

%At the two points beyond the boundary x values, h equals zero. 

H_EVAL = [0 H 0] ; 

%This is optimized code generated by MAPLE. SOL represents the 
%vector of solutions to this integral at each point x_i+l/2 and 
%x_i-l/2 multiplied by h. 

T1 = kw*H_EVAL; 

T2 = UW_AV G-U_AVG; 

T4 = H_EVAL."2 ; 

T5 = kf*T4; 
t12 = kf"2; 

T18 = l./(l + kw*H_EVAL/3 + tl2*T4/3 + kw*kf*T4.*H_EVAL/12); 

T20 = l./(l + Tl/3 + tl2*T4/3 + kw*kf*T4.*H_EVAL/12) ; 

T29 = T1.*T2.*T20; 

T35 = (1 + Tl/4)*kf.*T4.*U_AVG.*T20; 

T37 = -T29/6 - T35/6; 

T40 = Tl.*T2/2 - T5.*U_AVG/2 - Tl.*((l/3 + T5/12)*kw.*H_EVAL.* ... 

T2.*T20 + T5.*U_AVG.*T20/6)/2 - T5.*T37/2; 
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T41 = T4 0.~ 2; 

T45 = U_AVG - T29/6 - T35/6; 

T51 = T4.~2; 

T53 = T45.~ 2; 

T59 = T37."2; 

SOL = H_EVAL.*(T 41/5 + T5.*T45.*T40/2 +2/3*T37.*T40 + tl2*T51.* ... 
T53/3 + T37*kf.*T4.*T45 + T59); 

%The value of the integral over x at each point x is the difference 
%between adjacent terms 

PRES = .5*(([H 0]."2 - [0 H]."2) + diff(SOL)); 
else % resid = false case 

PRES = .5*(([H 0]."2 - [0 H ] . ~ 2 ) ) ; 

end 

%END PRESSURE TERMS 
POTENTIAL TERM 

BATH = feval(bath,X); %bathymetry evaluated at each point in X 
BATH_MID = BATH(1:end-1) + diff(BATH)/2; %bathymetry at midpoints in X 
POT = [0 H] .* (BATH - [0 BATH_MID]) + [H 0] .* ([BATH_MID 0] - BATH); 

%END POTENTIAL TERM 

FRICTION TERMS 

if friction == true 

%Evaluate wind at all X points from the user-defined function 
%passed to this method in wind. 

UW = feval(wind, X) ; 

%temporary solution: average h_i-l and h_i to get the value at x_i 
H_AVG = ([OH] + [H 0])/2; %This is may not be correct 

%This section calculates the values of u_tilde if the residual 
%functionality is turned on by the user, otherwise it returns a 
%zero array for the residual values, 
if resid == true 

%This section calculates the values of u_tilde at s=0 and s=l. 
%These formulae were automatically generated by Maple. 

T1 = kw*H_AVG; 

T2 = UW-U; 

T4 = H_AVG."2; 

T5 = kf*T4; 
t12 = kf"2; 

T18 = l./(l + kw*H_AVG/3 + tl2*T4/3 + kw*kf*T4.*H_AVG/12); 

T20 = l./(l + Tl/3 + tl2*T4/3 + kw*kf*T4.*H_AVG/12); 

UHI = (1/3 + T5/12)*kw.*H_AVG.*T2.*T18 + T5.*U.*T18/6; 

ULO = -Tl.*T2.*T20/6 - (1 + Tl/4)*kf.*T4.*U.*T20/6; 
else %resid = false case 

UHI = zeros(1, length(X)) ; 

ULO = UHI; 

end 


%This section creates the matrices M_fr and M_w used to multiply 
%the UW, U, U(l), and U(0) vectors. 

%Using this method means that extended edges can cause way too much 
%friction. 

%DELTA_X = diff(X)."3 

%By making the DELTA_X constant, we prevent having runaway friction 
%at the edges. 
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DELTA_X = ( (X(end)-X(l) ) / (n - 1))*3 * ones(l,n-l); 

MMWW = ((kw*mu)/(64*lambda"2))* spdiags([[DELTA_X 0]',... 

3*([0 DELTA_X] + [DELTA_X 0])',[0 DELTA_X]'],-1:1,n,n); 

DELTA_X_H = DELTA_X .* H; 

MMFF = ((kf*mu)/(64*lambda"2))* spdiags([[DELTA_X_H 0] ' ,... 

3*([0 DELTA_X_H] + [DELTA_X_H 0])', [ 0 DELTA_X_H]'], -1 : 1, n, n); 

%The simpler method is to pick a constand to represent M_fr and m_w 


end 

%END FRICTION 
solve for u 

%EXPLICIT, NO FRICTION 

%U = (MM\((MM + MMVV)*U' - dt*PRES' - dt*POT' )) ' ; 

if friction == false 

%IMPLICIT, NO FRICTION 

U = ((MM - MMVV)\(-dt*PRES' - dt*POT' + MM*U'))'; 

else 

%IMPLICIT, WITH FRICTION 

U = ((MM - MMVV + dt* (MMWW + MMFF)) \ (-dt*POT' - dt*PRES' + ... 
dt*MMWW*(UW - UHI)' - dt*MMFF*ULO' + MM*U'))'; 

end 

% update x at t_n+l 
X = X + U * dt; 

% calculate h at t_n+l 
H = M./diff(X); 

Plotting 

if i-T*floor(i/T) <1 | i == 1 


XX = X; 

HH = ([OH] + [H 0])/2; 

HH(1) = 0; 

HH(end) = 0; 

%calculate center of mass and momentum in continuing vectors 
%not a general solution... depends specifically on the masses being 
%equally distributed among the intervals 


XC = [XC sum(MM*X')/mass]; 

VC = [VC sum(MM*U')/mass]; 

%set title fields 

ti = ['Water Elevation - red. Bottom Profile - black. Water' ... 

' Height - blue vs. distance']; 
tstr = sprintf('Time: %.lf',t); 
istr = sprintf('Index: %d', i) ; 

indices = [tstr ' ' istr]; 

lstr = sprintf('Left Boundary: %.4f',X(l)); 

cmass = sprintf('Center of Mass: %.4f', XC(end)); 

rstr = sprintf('Right Boundary: %.4f',X(n)); 

dists = [Istr ' ' cmass ' ' rstr]; 
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%make second version of H for elevation above potential 
H_PLOT = HH + feval(bath, XX); 

%this is the lower end of the plot, this needs to be user set 

%somehow 

lower = 0; 


%Alternate displays for movie or otherwise 
if data.timing 

else if data.mov %movie case 
mov = mov + 1 

plot(X_POT,Y_POT,'k',XX,HH+lower,'b',XX,H_PLOT,'r',X(1), 
feval(bath,X(1)),'+ r', X(n),feval(bath,X (n)),'+ r' ) ; 
title(strvcat(ti, indices, dists)); 
axis([lo_b hi_b 0 ((3.5)~2)/2]) ; 

if i == 1 

set(1,'Position',[100,200,1000,400]) 
pause; 

playme(mov) = getframe; 


else 

playme(mov) = getframe; 

end 

else %not movie case 

%display height curves to first plot frame 
subplot(211); 

plot(X_POT,Y_POT,'k',XX,HH+lower,'b',XX,H_PLOT,'r',X(1), ... 

feval(bath,X(1)),'+ r',X(n),feval(bath,X(n)),'+ r'); 
title(strvcat(ti, indices, dists)); 
axis([lo_b hi_b 0 ((3.5)~2)/2] ) ; 

%axis([data.lo data.hi lower 6]); 

%display center of mass vs center of momentum plot 
subplot(212); 

plot(AX,AY,'k',AY,AX,'k',XC(2:end),VC(2:end),'b',XC(end) , 

VC(end),'+ b'); 

title('Center of Mass vs. Average Momentum - blue'); 
axis( [-2 2 -2 2]); 
axis square 

drawnow 
if data.mov 

playme(mov) = getframe; 

end 

if i == 1 
pause; 

end 

end 

end 


end 


end 

if data.timing 
toe 

end 


if data.mov 
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movie2avi(playme, data.name, 'fps', 10) 

end 


%%create the output data 
if data.timing 
out = toe; 

else 

out.res = resid; 

out.bath = bath; 

out.hi = hi_b; 

out.lo = lo_b; 

out.n = n; 

out. dt = dt ; 

out. X = X; 

out. H = H; 

out. U = lb- 

out, time = 0; 

out.k_hi = kw; 

out.k_lo = kf; 

out.uw = wind; 

out.lam = lambda; 

out.timing = data.timing; 

out.mov = data.mov; 

out.name = data.name; 

end 


end 
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D MATLAB Code: initialize3d.m 

function data = initialize3d() 


%% Data input for initial conditions 

% 2'“N+1 is the number of grid lines in each direction in the original 

% square domain 

N=5; 

% Coefficient of rotation: analogous to a coriolis force 
b=l.5; 

% Maximum time before computation stops 
Tmax=10; 

% Bottom friction coefficient: a=0 corresponds to no friction 
a = 0; 

% Initial x-coordinate of the fluid center of mass 
xc = 1/sqrt(2); 

% Initial y-coordinate of the fluid center of mass 
yc = 0; 

% Eddy viscosity 
mu= .05; 

% Ratio of dt/dx 
al = .01; 

%% Setting up domain and initial values of X, Y, u, and v 
m=2"N; 

x=-.5:1/m:.5; 
y=ones(size(1:m+l)); 

X=(y'*x) ' ; 

Y=(x'*y) ' ; 

dx=l/m ; 

H=1; 

dt=al*dx; 

d=0; 
s=0; 
w=0; 

T=0; 

u=((d+w)*X+(w+T)*Y)/2; 
v=((T-w)*X+(d-s)*Y)/2; 

X=X+xc; 

Y=Y+yc; 

%% Calculating area of the fluid elements 

DxXl(1:m,1:2:m-l)=X(2:m+l,1:2:m-l)-X(1:m,l:2:m-l); 

DxYl(l:m,1:2:m-1)=Y(2:m+1,1:2:m-1)-Y(1:m, l:2:m-l) ; 

DxXl(1:m,2:2:m)=X(2:m+1,3:2:m+1)-X(1:m,3:2:m+1); 

DxYl(1:m,2:2:m)=Y(2:m+l,3:2:m+1)-Y(1:m,3:2:m+1); 

DxX2 (1:m, 1:2 :m-1) =X (2 :m+1,2 : 2 :m) -X (1 :m, 2 : 2 :m) ; 

DxY2(l:m,1:2:m-1)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m); 
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DxX2(1: m, 2 :2 : m)=X(2:m+1, 2 : 2 : m)-X(1:m,2:2:m); 

DxY2(1:m,2:2:m)=Y(2:m+1, 2:2:m)-Y(1:m,2:2:m); 

DyXl(2:2:m,1:m)=X(2 : 2 :m, 2:m+l)-X(2:2:m,1:m); 

DyYl(2:2:m,1:m)=Y(2:2:m,2:m+l)-Y(2:2:m,1:m); 

DyXl(1:2:m-l,1:m)=X(2 : 2 :m, 2:m+l)-X(2:2:m, 1:m); 

DyYl(1:2:m-l,1:m)=Y(2 : 2 :m, 2:m+l)-Y(2:2:m,1:m); 

DyX2(1:2:m-l,1:m)=X(1:2:m-l,2:m+l)-X(1:2:m-l,1:m); 
DyY2(1:2:m-l,1:m)=Y(1:2:m-1,2:m+1)-Y(1:2:m-1,1:m); 
DyX2(2:2:m,1:m)=X(3:2:m+1,2:m+1)-X(3:2:m+1,1:m); 
DyY2(2:2:m,1:m)=Y(3:2:m+1,2:m+1)-Y(3:2:m+1,l:m); 


% area of type 1 and type 2 triangles assigned to each node 
Al=(DxXl.*DyYl-DyXl.*DxYl)/2; 

A2=(DxX2.*DyY2-DyX2.*DxY2)/2; 


% mass in each fluid element 

ml=H*Al; 

m2=H*A2; 


%% Mass Matrix Computation 

mass=sum(sum(sum(Al))+sum(sum(A2)))/m/m; 

M=zeros(m+1); 

M(1:m+1,1:m+1)=mass; 

% Each processor calculates its own mass matrix 
M(1,2:m)=mass/2; 

M(m+1,2:m)=mass/2; 

M(2:m,l)=mass/2; 

M(2:m,m+l)=mass/2; 

M(1,1)=mass/4; 

M(m+1,m+1)=mass/4; 

M(m+1,1)=mass/4; 

M(l,m+l)=mass/4; 

data.MM=sum(sum(M)) ; 

%% Constants for splitting to processors 
p = 3; 

P = 4; 

n = 2"(N-p+1); 


data.m = m; 
data.n = n; 
data.b = b; 
data.a = a; 
data.mu = mu; 
data.dt = dt; 
data.dx = dx; 
data.Tmax = Tmax; 
data.P = P; 
data.xc = xc; 
data.yc = yc; 

data.X = X; 
data.Y = Y; 
data.u = u; 
data.v = v; 
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data.M = M; 
data.ml = ml; 
data.m2 = m2; 

end 
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E MATLAB Code: shallowPar3d.m 


function plotData = shallow_par3d() 


%% Create a data input objects only on processor 1 
if labindex == 1 

input = initialize3d() ; 

P = input.P; 
m = input.m; 
n = input.n; 

pass(P).m = m; 

% constant data 
for i = 1:P 

pass(i).m = m; 

pass (i) .n = n+2; 

pass(i).b = input.b; 

pass(i).a = input.a; 

pass (i) .mu = input.mu; 

pass(i).dt = input.dt; 

pass(i).dx = input.dx; 

pass(i).runtime = input.Tmax; 

pass(i).P = input.P; 

pass(i).xc = input.xc; 

pass(i).yc = input.yc; 

end 

% vertex data 

pass(l).X = input.X(:,1:n+3); 
pass(l).Y = input.Y(:,1:n+3) ; 
pass(l).u = input.u(:, 1:n+3) ; 
pass(l).v = input.v(:,1:n+3); 
pass(l).M = input.M(:,1:n+3) ; 

for k = 2:input.P-1 

pass(k).X = input.X(:,(k-1)*n-l:k*n+l); 
pass(k) .Y = input.Y (:, (k-1)*n-l:k*n+l); 
pass(k) .u = input.u (:, (k-1)*n-l:k*n+l) ; 
pass(k) .v = input.v(:, (k-1)*n-l:k*n+l) ; 
pass(k).M = input.M(:,(k-1)*n-l:k*n+l); 

end 

pass(P).X = input.X(:,(P-1)*n-l:m+l); 
pass(P).Y = input.Y(:,(P-1)*n-l:m+l); 
pass(P).u = input.u(:,(P-l)*n-l:m+l); 
pass(P).v = input.v(:,(P-1)*n-l:m+1); 
pass(P).M = input.M(:,(P-1)*n-l:m+l); 

% triangle data 

pass(1).ml = input.ml(:,1:n+2); 
pass(1) .m2 = input.m2(:, 1:n+2) ; 

for k = 2:P-1 

pass (k) .ml = input.ml(:, (k-1)*n-l:k*n); 
pass(k) .m2 = input.m2(:, (k-1)*n-l:k*n) ; 

end 

pass(P).ml = input.ml(:,(P-1)*n-l:m); 
pass(P).m2 = input.m2(:,(P-1)*n-l:m); 
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end 


%% End setup, start distribution to labs 

if labindex > 1 

data = labReceive(1) ; 
else % in this case labindex == 1 
for i = 2:P 

labSend(pass(i) , i) ; 

end 

data = pass (1); 

end 


%% Set output object for constants 
if labindex == 1 

consts.M = input.M; 
consts.MM = input.MM; 
consts.xc = input.xc; 
consts.yc = input.yc; 
consts.m = input.m; 
consts.n = input.n; 

consts.plotpoints = input.Tmax*10; %this is not yet dynamic 

else 

consts = 0; 

end 


%% Begin run on each individual processor 

runtime = data.runtime; 
dt = data.dt; 

P = data.P; 

step=.l/dt; %used to determine when to plot 
k = 0; %index 

t = 0; %time elapsed 

%% Set up plotting variables 

plotData(runtime*10).empty = 0; %this is also not yet dynamic 

c = 0; %plot count 

%% Running Loop 
while t <= runtime 

k = k + 1; 

t = t + dt; 

%% Computation 

data = spOneStep(data); 

%% Communication with neighbors 

% Share data with adjacent processors 
n = data.n; 

if labindex == 1 

labSend([data.X(:,n-3:n-2) data.Y(:,n-3:n-2) data.u(:,n-3:n-2) ... 

data.v(:,n-3:n-2)] , 2) ; 
else if labindex == P 

fromLeft = labReceive(P - 1); 

else 

fromLeft = labSendReceive(labindex + 1, labindex - 1, ... 

[data.X(:,n-1:n) data.Y(:,n-1:n) data.u(:,n-1:n) ... 
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data.v(:,n-1:n)]); 

end 

end 

if labindex == 1 

fromRight = labReceive(labindex + 1); 
else if labindex == P 

labSend([data.X(:,3) data.Y(:,3) data.u(:,3) 

p-1) ; 

else if labindex == 2 

fromRight = labSendReceive(1, 3, ... 

[data.X(:,3:5) data.Y(:,3:5) data.u( 
data.v(:,3:5)]); 

else 

fromRight = labSendReceive(labindex - 1, 
[data.X(:,3) data.Y(:,3) data.u(:,3) 

end 

end 

end 

% Update matrices with data from adjacent processors 


if labindex 

'= i 



data.X( 

, i) 

= fromLeft ( 

,d 

data.Y( 

, i) 

= fromLeft( 

, 3 ) 

data.u( 

, i) 

= fromLeft( 

,5) 

data.v( 

, i) 

= fromLeft( 

, 1 ) 

data.X( 

, 2) 

= fromLeft( 

, 2 ) 

data.Y( 

, 2) 

= fromLeft( 

, 4 ) 

data.u( 

, 2) 

= fromLeft( 

, 6) 

data.v( 

, 2) 

= fromLeft( 

, 8 ) 


end 

if labindex ~= P 


if labindex 

== 1 




data.X( 

, n-1) 

= 

fromRight(:, 

1) ; 

data.Y( 

, n-1) 

= 

fromRight(:, 

4); 

data.u( 

, n-1) 

= 

fromRight(:, 

V) ; 

data.v( 

, n-1) 

= 

fromRight(:, 

10) 

data.X( 

, n) = 

fromRight(:,2) 

; 

data.Y( 

, n) = 

fromRight(:, 5) 

; 

data.u( 

, n) = 

fromRight(:,8) 

; 

data.v( 

, n) = 

fromRight(:,11); 

data.X( 

, n+1) 

= 

fromRight(:, 

3) ; 

data.Y( 

, n+1) 

= 

fromRight(:, 

6) ; 

data.u( 

, n+1) 

= 

fromRight(:, 

9) ; 

data.v( 

, n+1) 

= 

fromRight(:, 

12) 

else 





data.X( 

, n+1) 

= 

fromRight(:, 

1) 

data.Y( 

, n+1) 

= 

fromRight(:, 

2 ) 

data.u( 

, n+1) 

= 

fromRight(:, 

3 ) ; 

data.v( 

, n+1) 

= 

fromRight(:, 

4) ; 


end 

end 

Plotting every 1/10 second 
jj=k-step*floor(k/step); 
if jj == 0 

c = c + 1 


data.v(:,3)],... 

,3:5) ... 

labindex + 1, .. . 

data.v(:,3)]); 


if labindex 


1 
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plotData (c) .X 
plotData (c) .Y 
plotData(c).u 
plotData (c) .v 


data.X(:,1:end-3); 
data.Y(:,1:end-3); 
data.u(:,1:end-3) ; 
data.v(:,1:end-3) ; 


plotData(c) ,A1 = data,A1(:,1:end-3) ; 
plotData(c),A2 = data. A2(:,1:end-3); 


else if labindex == P 


plotData(c).X = 
plotData(c).Y = 
plotData(c).u = 
plotData(c).v = 


data.X(:, 3:end) ; 
data.Y(:,3:end) ; 
data.u(:, 3:end) ; 
data.v(:, 3:end) ; 


plotData(c) .A1 = data.A1(:, 2:end) ; 
plotData(c) ,A2 = data.A2(:,2:end) ; 


else 


plotData(c) .X = data.X(:,3:end-1) ; 
plotData(c) .Y = data.Y(:,3:end-1) ; 
plotData(c) .u = data.u(:,3:end-1) ; 
plotData(c).v = data.v(:,3:end-1); 


plotData(c) ,A1 = data.A1(:, 2:end-1) ; 
plotData(c) .A2 = data.A2(:,2:end-1) ; 


end 

end 

plotData(c).t = t; 


end 
%% Loop 

%% Assign objects into 'base' workspace for easy transfer to client 
assignin('base', 'out', plotData); 
assignin('base', 'static', consts) ; 


end 
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F MATLAB Code: spOneStep.m 

function data = spOneStep(data) 

%% Data in from last step 
% This data is constant 

m = data.m; % number of grid lines in x-direction 

n = data.n; % number of grid lines in y-direction 

b = data.b; % rotation coefficient 

a = data.a; % bottom friction coefficient 

mu = data.mu; % eddy viscocity 

dt = data.dt; % time step 

%dx = data.dx; 


ml = data.ml; % matrices of masses 
m2 = data.m2; 

M = data.M; 


% This data 
X = data.X; 
Y = data.Y; 
u = data.u; 
v = data.v; 


is changed 
% matrix of 
% matrix of 
% matrix of 
% matrix of 


x-coordinate 
y-coordinate 
x-component s 
y-component s 


values 
values 
of velocity 
of velocity 


%% Calculate properties along edges of triangles 
% positions differences 

DxXl(1 :m,1:2:n-1)=X(2:m+l,1:2:n-1)-X(1:m,l:2:n-l); 
DxYl(l:m,1:2:n-1)=Y(2:m+l,1:2:n-1)-Y(1:m,l:2:n-l) ; 
DxXl(1:m,2:2:n)=X(2:m+1,3:2:n+1)-X(1:m,3:2:n+1); 
DxYl(1:m,2:2:n)=Y(2:m+l,3:2:n+1)-Y(1:m,3:2:n+1); 

DxX2(1: m, 1:2 :n-1)=X(2:m+1,2:2:n)-X(1 :m, 2 : 2 : n); 

DxY2(1:m,1:2:n-l)=Y(2:m+l,2:2:n)-Y(1:m, 2:2:n) ; 

DxX2(1:m,2:2:n)=X(2:m+1,2:2:n)-X(1:m,2:2:n); 

DxY2(1:m,2:2:n)=Y(2:m+1, 2:2:n)-Y(1:m,2:2:n); 

DyXl(2:2:m,1:n)=X(2:2:m,2:n+1)-X(2:2:m,1:n); 

DyYl(2:2:m,1:n)=Y(2:2:m,2:n+l)-Y(2:2:m,1:n); 

DyXl(1:2:m-1,1:n)=X(2:2:m,2:n+1)-X(2:2:m,1:n); 

DyYl(1:2:m-l,1:n)=Y(2:2:m,2:n+1)-Y(2 : 2 :m, 1: n) ; 

DyX2(1:2:m-l,1:n)=X(1:2:m-l,2:n+1)-X(1:2:m-l,1:n); 
DyY2(1:2:m-l,1:n)=Y(1:2:m-l,2:n+1)-Y(1:2:m-l,1:n); 

DyX2(2:2:m,1:n)=X(3:2:m+1,2:n+1)-X(3:2:m+1, 1:n) ; 
DyY2(2:2 :m,1:n)=Y(3:2:m+l,2:n+1)-Y(3:2:m+1,l:n); 

% calculate this data to avoid storing it 

Al=(DxXl.*DyYl-DyXl.*DxYl)/2; 

A2=(DxX2.*DyY2-DyX2.*DxY2)/2; 

hl=ml./Al; 
h2=m2./A2; 

% velocity differences 

Dxul(1:m,1:2:n-1)=u(2:m+1,1:2:n-1)-u(1:m, 1:2:n-1) ; 
Dxvl(l:m,1:2:n-1)=v(2:m+l,1:2:n-1)-v(1:m,l:2:n-l); 
Dxul(1:m,2:2:n)=u(2:m+1,3:2:n+1)-u(1:m,3:2:n+1); 
Dxvl(1:m,2:2:n)=v(2:m+1,3:2:n+1)-v(1:m,3:2:n+1); 
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Dxu2 (1:m,1:2:n-1) =u (2 :m+l ,2:2:n)-u(l:m,2:2:n) ; 

Dxv2(l:m,l:2:n-l)=v(2:m+l, 2:2 : n)-v(1:m, 2 :2 : n); 

Dxu2 (1: m, 2 :2 : n) =u (2 : m+1,2:2:n)-u(l:m,2:2:n) ; 

Dxv2(1:m,2:2:n)=v(2:m+1, 2:2 : n) -v (1: m, 2 :2 : n); 

Dyul (2 : 2 :m, 1: n) =u (2 : 2 :m, 2 :n+1) -u (2 : 2 :m, 1:n) ; 

Dyvl (2 : 2 :m, 1: n) =v (2 : 2 :m, 2 :n+1) -v (2 : 2 :m, 1:n) ; 

Dyul(1:2:m-l,1:n)=u(2:2:m,2:n+l)-u(2:2:m, l:n) ; 

Dyvl(1:2:m-l,1:n)=v(2:2:m,2:n+l)-v(2:2:m,1:n); 

Dyu2(1:2:m-l,1:n)=u(1:2:m-1,2:n+1)-u(1:2:m-l, 1:n) ; 

Dyv2(1:2:m-1,1:n)=v(1:2:m-1,2:n+1)-v(1:2:m-1,1:n); 

Dyu2 (2 : 2 : m, 1: n) =u (3:2 :m+l, 2 : n+1) -u (3 : 2 : m+1, 1: n) ; 

Dyv2(2:2:m,1:n)=v(3:2:m+l,2:n+l)-v(3:2:m+l, 1:n) ; 

%% Pressure terms in x momentum eqn. 

plRLx=(hi. ~ 2 ). *(DyYl)/4; 
p2RLx=(h2. ~ 2 ). *(DyY2)/4; 
plUDx=(hi. ~ 2 ). *(DxYl)/4; 
p2UDx=(h2."2).*(DxY2)/4; 

Plrxx=0*X; 

Plrxx(1:m,1:n)=plRLx; 

PlLxx=0*X; 

PILxx(2:m+1 , 1:n)=plRLx; 

P2rxx=0*X; 

P2rxx(1:m,1:n)=p2RLx; 

P2Lxx=0*X; 

P2Lxx(2:m+1 , 1:n)=p2RLx; 

P2Uxy=0*Y; 

P2Uxy(1:m,1:n)=p2UDx; 

P2Dxy=0*Y; 

P2Dxy(1:m,2:n+1)=p2UDx; 

PlUxy=0*Y; 

PlUxy(1:m,1:n)=plUDx; 

PlDxy=0*Y; 

PIDxy(1:m,2:n+1) = plUDx; 

% Pressure Forces in x eqn 
PPx=0*X; 

PPx(:,1)=PPx(:,1)+ dt*(PILxx(:,1)-Plrxx(:,1)); 

PPx(:,3:2:n+1)=PPx(:,3:2:n+1)+dt*( PILxx(:,3:2:n+1)-Plrxx(:,3:2: 
PPx(:,3:2:n+1)=PPx(:,3:2:n+1)+dt*(PILxx(:,2:2:n)-Plrxx(:,2:2:n)) 
PPx(:,2:2:n)=PPx(:,2:2:n)+dt*(P2Lxx(:,1:2:n-1)-P2rxx(:,1:2:n-1)) 
PPx(:,2:2:n)=PPx(:,2:2:n)+dt*(P2Lxx(:,2:2:n)-P2rxx(:, 2:2:n)) ; 

PPx(1,:)=PPx(1,:)-dt*(P2Dxy(1,:)-P2Uxy(1,:)); 

PPx(3:2:m+1,:)=PPx(3:2:m+l,:)-dt*(P2Dxy(3:2:m+l,:)-P2Uxy(3:2:m+l 
PPx(3:2:m+1,:)=PPx(3:2:m+l,:)-dt*(P2Dxy(2:2:m,:)-P2Uxy(2:2:m,:)) 
PPx(2:2:m,:)=PPx(2:2:m,:)-dt*(PIDxy(1:2:m-l,:)-PlUxy(1;2:m-l,:)) 
PPx(2:2:m,:)=PPx(2:2:m,:)-dt*(PIDxy(2:2:m,:)-PlUxy(2:2:m,:)); 


n+1) ) ; 




%% Pressure terms in y momentum eqn. 



p2udy=(h2. ~ 2 ). *(DxX2)/4; 
pludy=(hi. ~ 2 ). *(DxXl)/4; 
p2RLy=(h2. ~ 2 ). *(DyX2)/4; 
plRLy=(hi. ~ 2 ). *(DyXl)/4; 

P2uyy=0*Y; 

P2uyy(1:m,1:n)= p2udy; 
P2dyy=0*Y; 

P2dyy(1:m,2:n+1)=p2udy; 
Pluyy=0*Y; 

Pluyy(1:m,1:n)= pludy; 
Pldyy=0*Y; 

Pldyy(1:m,2:n+l)= pludy; 


Plryx=0*X; 

Plryx(1:m,1:n)= plRLy ; 
PlLyx=0*X; 

PILyx(2:m+l,1;n)=plRLy; 
P2ryx=0*X; 

P2ryx(1:m,1;n)=p2RLy; 
P2Lyx=0*X; 

P2Lyx(2:m+1,1:n)= p2RLy; 


% Pressure Forces in y eqn 
PPy=0*Y; 

PPy(1,:)=PPy(1,:)+dt*(P2dyy(1,:)-P2uyy(1,:)); 

PPy(3:2:m+1, :)=PPy(3:2:m+1,:)+dt*(P2dyy(3:2:m+l,:)-P2uyy(3:2:m+l,:)); 
PPy(3:2:m+1,:)=PPy(3:2:m+l,:)+dt*(P2dyy(2:2:m,:)-P2uyy(2:2:m,:)); 

PPy(2:2:m,:)=PPy(2:2:m,:)+dt*(Pldyy(1:2:m-l,:)-Pluyy(1:2:m-1,:)); 

PPy(2:2:m,:)=PPy(2:2:m,:)+dt*(Pldyy(2:2:m,:)-Pluyy(2:2:m,:)); 

PPy(:,1)=PPy(:,1)-dt*(PILyx(:,1)-Plryx(:,1)); 

PPy(:,3:2:n+1)=PPy(:,3:2:n+1)-dt*(PILyx(;,3:2:n+1)-Plryx(:,3:2:n+l)); 
PPy(:,3:2:n+1)=PPy(:,3:2:n+1)-dt*(PILyx(;,2:2:n)-Plryx(:,2:2:n)); 

PPy(:,2:2:n)=PPy(:,2:2:n)-dt*(P2Lyx(:,2:2:n)-P2ryx(:,2:2:n)); 

PPy(:, 2:2:n)=PPy(:,2:2:n)-dt*(P2Lyx(:,1:2:n-1)-P2ryx(:,1:2:n-1)); 


% this ends the pressure computation 
%% viscous stresses computation 

sl=mu*(hi).*(-(DxXl).*(Dyvl)+(Dxul).*(DyYl)+(DyXl).*(Dxvl)-(DxYl).*(Dyul)) 
s2=mu*(h2).*(-(DxX2).*(Dyv2)+(Dxu2).*(DyY2)+(DyX2).*(Dxv2)-(DxY2).*(Dyu2)) 

Tl=mu*(hi).*((DxXl).*(Dyul)+(DyYl).*(Dxvl)-(DyXl).*(Dxul)-(DxYl).*(Dyvl)); 
T2=mu*(h2).*((DxX2).*(Dyu2)+(DyY2).*(Dxv2)-(DyX2).*(Dxu2)-(DxY2).*(Dyv2)); 


Vlrxx=0*X; 

Vlrxx(1:m,1:n)=(si.*DyYl-Tl.*DyXl)/2; 



VlLxx=0*X; 

VILxx(2:m+l,l:n)=(si.*DyYl-Tl.*DyXl)/2; 


V2rxx=0*X; 

V2rxx(1:m,1:n)=(s2.*DyY2-T2.*DyX2)/2; 
V2Lxx=0*X; 

V2LXX(2:m+l,l:n)=(s2.*DyY2-T2.*DyX2)/2 ; 


Vluxy=0*X; 

Vluxy(l:m,1:n)=(T1.*DxXl-sl.*DxYl)/2; 
Vldxy=0*X; 

Vldxy(1:m,2:n+l)=(T1.*DxXl-sl.*DxYl)/2; 
V2uxy=0*X; 

V2uxy(l:m,1:n)=(T2.*DxX2-s2.*DxY2)/2; 

V2dxy=0*X; 

V2dxy(1:m,2:n+l)=(T2.*DxX2-s2.*DxY2)/2; 


% Viscous Forces in x eqn 
Wx=0*X; 

Wx ( : , 1) =Wx ( :, 1) + (Vlrxx ( : , 1) -VILxx ( : , 1) ) ; 

Wx (:, 3 : 2 : n+1) =VVx (:,3:2:n+l) + (Vlrxx (:,3:2:n+l) -VILxx (:, 3:2 :n+l) ) 
Wx {:, 3 : 2 : n+1) =VVx (:, 3 :2 : n+1) + (Vlrxx (:, 2:2 : n) -VILxx (:, 2 :2 : n) ) ; 

Wx (:, 2 : 2 : n) =Wx (:, 2 : 2 : n) + (V2rxx (:, 2 : 2 : n) -V2Lxx (:, 2 : 2 : n) ) ; 

Wx (:, 2 : 2 : n) =Wx (:, 2 : 2 : n) + (V2rxx (:, 1: 2 : n-1) -V2Lxx (:, 1:2 : n-1) ) ; 

Wx (1, : ) =VVx (1, : ) + (V2uxy (1, : ) -V2dxy (1, : ) ) ; 

Wx (2 : 2 :m, : ) =Wx (2 : 2 :m, : ) + (Vluxy (2 : 2 :m, : ) -Vldxy (2:2:m, :) ) ; 

Wx (2 : 2 :m, :) =Wx (2 : 2 :m, :) + (Vluxy (1: 2 :m-l, :) -Vldxy (1: 2 :m— 1, :) ) ; 

Wx (3:2:m+1, : ) =VVx (3:2:m+1, : ) + (V2uxy (3:2 :m+l, : ) -V2dxy (3:2 :m+l , : ) ) 
Wx (3:2:m+1 , :) =VVx (3:2:m+1 , :) + (V2uxy (2 :2 :m, :) -V2dxy (2 :2 :m, : ) ) ; 


Vlryx=0*X; 

Vlryx(1:m,1:n)=(T1.*DyYl+sl.*DyXl)/2; 
VlLyx=0*X; 

VILyx(2:m+l,1:n)=(T1.*DyYl+sl.*DyXl)/2; 
V2ryx=0*X; 

V2ryx(1:m,1:n)=(T2.*DyY2+s2.*DyX2)/2; 
V2Lyx=0*X; 

V2Lyx(2:m+l,1:n)=(T2.*DyY2+s2.*DyX2)/2; 


Vluyy=0*X; 

Vluyy(1:m,1:n)=(T1.*DxYl+sl.*DxXl)/2; 


Vldyy=0*X; 
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Vldyy(1:m,2:n+l)=(T1.*DxYl+sl.*DxXl)/2; 
V2uyy=0*Y; 

V2uyy(1:m,1:n)=(T2.*DxY2+s2.*DxX2)/2; 
V2dyy=0*X; 

V2dyy(1:m,2:n+l)=(T2.*DxY2+s2.*DxX2)/2; 


% Viscous Forces in y eqn 
Wy=0*Y; 

Wy (:, 1) =Wy (:, 1) + (Vlryx ( :, 1) -VILyx (:, 1) ) ; 

Wy (: , 3 : 2 : n+1) =VVy (: , 3 : 2 : n+1) + (Vlryx (: , 3 :2 : n+1) -VILyx (:, 3 : 2 : n+1) ) ; 
Wy (:, 3 : 2 : n+1) =VVy (:, 3 : 2 : n+1) + (Vlryx (:, 2 : 2 : n) -VILyx (:, 2 :2 : n) ) ; 

Wy (:, 2 : 2 : n) =Wy (:, 2 : 2 : n) + (V2ryx (:, 2 : 2 : n) -V2Lyx (:, 2 : 2 : n) ) ; 

Wy (:, 2 : 2 : n) =Wy (:, 2 : 2 : n) + (V2ryx (:, 1: 2 : n-1) -V2Lyx (:, 1: 2 : n-1) ) ; 

Wy (1, :) =Wy (1, :) - (V2uyy (1, :) -V2dyy (1, :) ) ; 

Wy (3:2 :m+l, :) =VVy (3:2 :m+l, :) - (V2uyy (3:2 :m+l, :) -V2dyy (3:2 :m+l, :) ) ; 
Wy (3:2 :m+l, :) =VVy (3:2 :m+l, :) - (V2uyy (2 : 2 :m, :) -V2dyy (2 : 2 :m / :) ) ; 

Wy (2 : 2 :m, :) =Wy (2 :2 :m, :) - (Vluyy (2 :2 :m, :) -Vldyy (2 :2 :m, :) ) ; 

Wy (2 : 2 :m, :) =Wy (2 :2 :m, :) - (Vluyy (1: 2 :m-l, :) -Vldyy (1: 2 :m-l, :) ) ; 

% this concludes viscous stress computation 

%% Computation of x and y Forces 

Fx=PPx+Wx; 

Fy=PPy+Wy; 


% this concludes the force computation 

%% Velocity update 

u=(l-a*dt)*u+(Fx./M-dt*(X-b*v)) ; 
v=(l-a*dt)*v+(Fy./M-dt*(Y+b*u)) ; 


%% Leap Frog 

X=X+dt*u; 
Y=Y+dt*v; 


%% Store data to pass to next step 

%update some data 

Al= (DxXl.*DyYl-DyXl.*DxYl)/2; 

A2=(DxX2.*DyY2-DyX2.*DxY2)/2; 


% only 

store 

data that 

data.X = 

X; 

% 

matrix 

of 

data.Y = 

Y; 

% 

matrix 

of 

data.u = 

u; 

% 

matrix 

of 

data.v = 

v; 

% 

matrix 

of 


has been changed 
x-coordinate values 
y-coordinate values 
x-components of velocity 
y-components of velocity 


data.Al = Al; 
data.A2 = A2; 
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% This data is stored for plotting - DEPRECATED FOR NOW 
%SS = si.*sl+Tl.*Tl+s2.~2 +T2."2; 

%data.SSl_l = max(max(SS)); 

%data.SSl_2 = max(max(max(A1)),max(max(A2)))/dx/dx; 
%data.SSl_3 = min(min(A2))/dx/dx; 

end 
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G MATLAB Code: spPlot.m 

function spPlotO 

%% Load data from labs to client 
pmode lab2client out 1 setl 
pmode lab2client out 2 set2 
pmode lab2client out 3 set3 
pmode lab2client out 4 set4 
pmode lab2client static 1 static 

%% Move data from base workspace to function workspace 

static = evalin('base', 'static'); 

setl = evalin('base', 'setl'); 

set2 = evalin('base', 'set2'); 

set3 = evalin('base', 'set3'); 

set4 = evalin('base', 'set4'); 

%% Extract data from structures 
xc = static.xc; 
yc = static.yc; 

Time = 0; 

R = xc"2+yc"2; 

m = static.m; 
n = static.m; 

M = static.M; 

MM = static.MM; 

plotpoints = static.plotpoints; 

[pData(1:plotpoints).empty] = deal(0); 
for i = liplotpoints 

pData(i).X = [setl(i).X set2(i).X set3(i).X set4(i).X]; 
pData(i).Y = [setl(i).Y set2(i).Y set3(i).Y set4(i).Y]; 
pData(i).u = [setl(i).u set2(i).u set3(i).u set4(i).u]; 
pData(i).v = [setl(i).v set2(i).v set3(i).v set4(i).v]; 

pData(i).Al = [setl(i).Al set2(i).Al set3(i).Al set4(i).Al]; 
pData(i).A2 = [setl(i).A2 set2 (i) .A2 set3(i).A2 set4(i).A2]; 

pData(i).t = setl(i).t; 

end 

%% Plot 

for c = l:plotpoints 

X = pData(c).X; % matrix of x-coordinate values 

Y = pData(c).Y; % matrix of y-coordinate values 

u = pData(c).u; % matrix of x-components of velocity 

v = pData(c).v; % matrix of y-components of velocity 

A1 = pData(c). Al; 

A2 = pData(c).A2; 

t = pData(c).t; 

HH=0*X; 

AA=0*X; 

AA(2:2:m,2:2:n)=(A2(2:2:m,2:2:n)+A2(2:2:m,1:2:n-l)+A2(1:2:m-l,2:2:n)+... 

A2(1:2:m-l,1:2:n-l))/4; 



AA (2:2:m,2:2:n) =AA(2 :2 :m, 2 : 2 : n) + (A1 (2 : 2 : m, 2 : 2 : n) +A1 (2 :2 :m, 1:2:n-l) + . . . 

A1(1:2:m-1,2:2:n)+...A1(1:2:m-l,1:2:n-l))/4 
AA(2 : 2 :m, 3 :2 : n-1) = (A1(2:2:m,3:2:n-l)+A1(1:2:m-l,3:2:n-1)+... 

A1(1:2:m-l,2:2:n-2)+A1(2:2:m,2:2:n-2))/2; 

AA(3:2:m-l,2:2:n)=(A2(3:2:m-l,2:2:n)+A2(3:2:m-l,l:2:n-l)+... 

A2(2:2:m-2,2:2:n)+A2(2:2:m-2,l:2:n-l))/2; 

AA(3:2:m-l,3:2:n-1)=(A2(3:2:m-l,3:2:n-1)+A2(2:2:m-2,3:2:n-1)+... 

A2(3:2:m-l,2:2:n-2)+A2(2:2:m-2,2:2:n-2))/4; 
AA(3:2:m-l,3:2:n-1)=AA(3:2:m-1,3:2:n-1)+(A1(3:2:m-l,3:2:n-1)+... 

A1(1:2:m-2,3:2:n-l)+A1(3:2:m-l,2:2:n-2)+A1(2:2:m-2,2:2:n-2))/4 
HH (2 :m, 2 : n) =M (2 :m, 2 : n) ./AA(2:m,2:n) ; 


Xc=sum(sum(M. *X))/MM; 

Yc=sum(sum(M.*Y))/MM; 
uc=sum(sum(M.*u))/MM; 
vc=sum(sum(M.*v))/MM; 

% This will do for now 
xc=[xc,Xc]; 
yc=[yc,Yc]; 

if max(size(R))<100; 

R=[R,Xc~2+Yc~2]; 

Time=[Time,t]; 

else 

R=[R(2:100),Xc~2+Yc~2]; 
Time=[Time(2:100),t]; 

end 


XX=X; 
YY=Y; 


ff=('velocity fields u-uc and v-vc — both approximately planar fields'); 


ffl=['surface plots 
cl f 


of water elevation h+(x"2+y"2)/2 
boundary, and time ' 


the parabolic ... 
num2str(t)] ; 


subplot(2,2,1 ) 
surf(XX,YY,HH) 
shading interp 
hold on 

contour3(XX,YY,HH,20,'k') 
hold on 

axis([-22-2202]) 
axis square 

title('surface plot of water column above z=(x"2+y"2)/2') 

subplot(322) 
surf(XX-Xc,YY-Yc,u-uc) 
shading interp 
hold on 

surf(XX-Xc,YY-Yc,v-vc) 
shading interp 

contour3(XX-Xc,YY-Yc,u-uc, 20, ' k' ) 
title (ff) 

axis([-1.5 1.5 -1.5 1.5 -1.5 1.5])% -.5 .5]) 
axis square 


subplot(3,2,4) 

plot(xc,yc,'r',Xc,Yc,'+k',0,0,'+k',[-2 2],[0 0],'k',[0 0],[-2 2],'k') 
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axis([-1 1-11]) 

title ('trajectory of the center of mass-initial position (l/sqrt(2) , 0) ') 

axis square 

subplot(223) 

ZZ1=((XX.~2+YY."2)/2) ; 

ZZ2=HH+(XX."2+YY.~2)/2; 
surf(XX,YY,ZZ1) 
shading interp 
hold on 

surf(XX,YY,ZZ2) 
shading interp 
hold on 

contour3(XX,YY,ZZ2,20,'k') 
hold on 

contour(XX,YY,ZZ2,20,'k') 
axis([-22-2202.5]) 
axis square 
title (ffl) 

subplot(3,2,6) 
plot (Time, R) 

axis([min(Time),max(Time), min(R),max(R)]) 

title (' xc"2+yc / '2 vs time') 

drawnow 

%pause(0.1) 


end 



